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Semi-analytic reconstruction uses the analytic solution to a second-order, steady, 
ordinary differential equation (ODE) to simultaneously evaluate the convective and 
diffusive flux at all interfaces of a finite volume formulation. The second-order ODE 
is itself a linearized approximation to the governing first- and second- order partial 
differential equation conservation laws. Thus, semi-analytic reconstruction defines 
a family of formulations for finite volume interface fluxes using analytic solutions to 
approximating equations. Limiters are not applied in a conventional sense; rather, 
diffusivity is adjusted in the vicinity of changes in sign of eigenvalues in order to 
achieve a sufficiently small cell Reynolds number in the analytic formulation across 
critical points. Several approaches for application of semi-analytic reconstruction 
for the solution of one-dimensional scalar equations are introduced. Results are 
compared with exact analytic solutions to Burger’s Equation as well as a conven- 
tional, upwind discretization using Roe’s method. One approach, the end-point 
wave speed (EPWS) approximation, is further developed for more complex appli- 
cations. One-dimensional vector equations are tested on a quasi one-dimensional 
nozzle application. The EPWS algorithm has a more compact difference sten- 
cil than Roe’s algorithm but reconstruction time is approximately a factor of four 
larger than for Roe. Though both are second-order accurate schemes. Roe’s method 
approaches a grid converged solution with fewer grid points. Reconstruction of flux 
in the context of multi-dimensional, vector conservation laws including effects of 
thermochemical nonequilibrium in the Navier- Stokes equations is developed. 


I. Nomenclature 

Bold face, lowercase variable names refer to vectors. Bold face, uppercase variable names refer 
to matrices. 


Roman symbols 

a 


A 

A 



c 


Cs 

Cp 

C, 

Ds 


constant in exact solution to Burgers Eq. used in Eqs. 51 and 53 
area 

inviscid-flux Jacobian with respect to conserved variables, df/dq 
viscous-flux Jacobian with respect to gradient of primitive variables, dhi/d-^^ 

speed of sound, + P)p/ p 
mass fraction of species s 
heat capacity at constant pressure 
heat capacity at constant volume 
mass diffusion coefficient for species s 
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Roman symbols, 

e 

E 

f,/ 
h, h 

hs,tr 

hs^v 

H 

k 

L2 

Ixjlyj Iz 

mx,my,mz 

M 

M 

Ms 

1 fly j IT-z 

P 

Pr 

Ppr 

Pr^ 

q,9 

q',Q' 

q 

ffMO, ffMl 

R 

Re 

Rcl 

Rcr 

Rex 

s, s 

Sc, 

t 

T 

T 

U, V, w 

U,V,W 

x,y,z 

X,Y,Z 


continued 

energy per unit mass 

total energy per unit mass, e + d 

inviscid flux (vector, scalar) 

viscous flux (vector, scalar) 

translational + rotational enthalpy for species s 

vibrational + electronic enthalpy for species s 

total enthalpy, E + p/p 

thermal conductivity 

error norm, sum of squares 

components of unit vector in Y direction 

components of unit vector in Z direction 

Mach number 

mixture molecular weight 

molecular weight of species s 

components of unit vector in X direction 

pressure 

Prandtl number, pCp/fj 

Prandtl number for translational + rotational energy, pCp/fjtr 

Prandtl number for vibrational + electronic energy, pCp/fjv 

conserved variable (vector, scalar) 

characteristic variable (vector, scalar) 

primitive variable vector 

reconstruction factors, Eq. 40 

matrix of right eigenvectors of A 

reference cell Reynolds number for edge, Eq. 45 

cell Reynolds number referenced to node L, Eq. 39 

cell Reynolds number referenced to node R, Eq. 39 

cell Reynolds number, 

source term (vector, scalar) 

Schmidt number for species s, ppDg 
time 

temperature 

transform from primitive to conservative variables, Tdq = dq 

Cartesian velocity components 

velocity components in rotated system 

reference Cartesian coordinates 

rotated Cartesian system 
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Greek symbols 

a, f3 
a 

P 

Is 

V 

V 

6 

A, A 

V 

P 

Ps _ 

4 >, 4 > 

$ 

Vi 

r 

tjj 

Subscripts 

1,2,3 referring to X, X, Z directions, respectively 

A, B, C, D nodes of tetrahedron 

CW S constant wave speed formulation 

e exact 

i,j node index or vector element 

L left node 

M dummy index for left or right node 

R right node 

Roe Roe formulation 

tr translational + rotational energy component 

V vibrational + electronic energy component 

II. Introduction 

The evolution of structured grid to unstructured-grid techniques for flow simulation is motivated 
by the relative ease of generating an initial grid over a complex configuration and the relative 
flexibility in adapting to flow features.^ Feature adaptation is particularly important in hypersonic 
flow simulation. The interaction of shocks, boundary layers, free shear-layers, and vortices spawn 
ever more complex flow structures at oblique angles to flow boundaries and upstream feature 
orientation. Such interactions are encountered in scram-jet engines, across compression surfaces, 
protuberances, and deflected control surfaces, around reaction control system jets, and even in 
the wake of simple blunt bodies. Inadequate resolution of these structures in a flow simulation 
adversely impacts the prediction of aerodynamic forces and surface heating. 


transformed coordinates referenced to left and right states, Eq. 30 
kinetic energy per unit mass, + v'^ + w'^) 
a perfect gas) 

dp 

dps 

transformed coordinate, Eq. 7 
thermal conductivity 

1 if mass fraction driven binary diffusion, M/Mg if mole fraction driven 

eigenvalue of A (diagonal matrix, scalar) 

viscosity 

diffusivity, /j,/p 

density 

density of species s 

artificial diffusivity factor, see Eqs. 46 - 48 
diagonal matrix operating on s to get analytic solution, Eq. 63 
diagonal matrix operating on dq/dx to get analytic solution, Eq. 62 
element of T 

diagonal part of R~^BT“^R 

extrapolation factor defining neighborhood of eigenvalue sign change used in Eq. 48 
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However, before unstructured simulations can be routinely applied to hypersonic flow simula- 
tions, a serious deficiency in simulation quality must be corrected. Unstructured grid simulations of 
hypersonic flows poorly predict stagnation region heating using high aspect ratio tetrahedra across 
the boundary layer. The deficiency is readily apparent in a simple test problem of hypersonic flow 
over a cylinder in which a biased, unstructured grid is used that includes spanwise resolution.^ In 
essence, the grid admits 3D flow for a case in which 2D flow is expected. Some jags in the shock 
capture process induce irregular entropy distribution in streamlines heading toward the stagnation 
point. On top of this non- uniformity, the element diagonals obliquely cut across the boundary layer 
inducing cross flow. The simulated flow quality is highly dependent on the flow limiters applied 
within the algorithm but, so far, there has been no algorithm fix identified which cures the problem 
in the context of the underlying, quasi one-dimensional reconstruction technique. The problem 
can be masked through use of prisms across the boundary layer with faces parallel or orthogonal 
to the gradient vector. However, unless one is prepared to abandon tetrahedra and automatically 
introduce and align prisms in all free shear layers and wakes then one must seek an algorithmic 
solution. 

The one-dimensional nature of reconstruction in current finite-volume methods is suspected to 
be a major contributor to the problem cited above. Reconstruction of flux across a cell face is based 
on a Riemann problem for a ID discontinuity cross the cell face.^^® The reconstruction process does 
not consider the principle flow direction - it is completely a function of the local grid orientation. 
In this sense, a structured grid has advantages over an unstructured grid crossing a single shock 
or boundary layer. All surfaces of the structured grid can be made parallel or orthogonal to the 
principle flow direction and gradient vectors. A tetrahedron will always have at least one edge 
oblique to the flow direction and gradient vector but is constrained to use this direction in the flux 
reconstruction process. 

This inherent difficulty with tetrahedra was a primary motivation to develop an algorithm 
for multidimensional flux reconstruction. The semi-analytic reconstruction introduced here uses 
a solution to the linearized conservation equations to define the variation of dependent variables 
between cells. Application of a semi- analytic approach is also inspired by success achieved with 
a post-processing algorithm that couples analytic solutions to CFD using an integral boundary 
layer approach for analyzing global changes to surface catalysis on heating^ and on a triple deck 
approach®’ ® for analyzing local changes to surface boundary conditions on heating. The challenge to 
the semi-analytic solution approach is the assumption of a continuous solution, even in the presence 
of shocks, that implicitly requires physical dissipation be included in the derivation. Effects of 
source terms on the reconstructed ffux are also included in the derivation; an effect usually ignored 
in finite- volume reconstruction although it has been discussed previously by Van Leer.^® 

Semi-analytic reconstruction developed herein uses the analytic solution to a second-order, 
steady, ordinary differential equation (ODE) to simultaneously evaluate the convective and diffu- 
sive ffux at all interfaces of a finite volume formulation. The second-order ODE is itself a linearized 
approximation to the governing first- and second- order partial differential equation (PDE) conser- 
vation laws. Thus, semi-analytic reconstruction defines a family of formulations for finite volume 
interface fluxes using analytic solutions to approximating equations. Limiters are not applied in a 
conventional sense; rather, diffusivity is adjusted in the vicinity of changes in sign of eigenvalues 
in order to achieve a sufficiently small cell Reynolds number in the analytic formulation across 
critical points. Several approaches for application of semi- analytic reconstruction for the solution 
of one-dimensional scalar equations are introduced in Sec. HI. Results are compared with exact 
analytic solutions to Burger’s Equation as well as a conventional, upwind discretizations using 
Roe’s scheme. One approach, the end-point wave speed (EPWS) approximation, is further devel- 
oped in Sec. IV for more complex applications. One-dimensional vector equations are tested on a 
quasi one-dimensional nozzle application. These test problems include strong shocks and source 
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terms. Possible approaches to higher-order accurate approximations are suggested in the Appendix, 
though these ideas have not been tested on anything more than a ID scalar equation. 

The ultimate proof of the algorithm in 3D will be the cylinder problem discussed above. These 
tests have not been completed but reconstruction of flux in the context of multi-dimensional, vector 
conservation laws including effects of thermochemical nonequilibrium in the Navier-Stokes equations 
is developed in Sec. V for future testing. 


III. One-Dimensional Scalar Equation 


A one-dimensional, scalar conservation law including both convection and diffusion can be 
expressed by the following equation. 


da d 


(1) 


A numerical solution to Eq. 1 written in conservative form must guarantee that the convective 
(/)and diffusive {h) flux exiting cell i is exactly balanced by that entering cell i + 1. Such a balance 
is enforced in the following difference equation. 


„n+l 


At 


+ 


(/ ^)i+l/2 (/ ^)i-l/2 


{Xi+l /2 — Xi-112) 


= Si 


( 2 ) 


Eq. 2 requires a reconstruction of / and h on cell boundaries. In a typical formulation the convective 
term and diffusive terms are reconstructed independently. Eor example, in Roe’s method, the 
convective term is defined as follows. 


/i+l/2 — 2 /*+! + /* “ l-^i-l-1/21 (qi+I - Qi + dq^J^/2) 


( 3 ) 


The variable Aj_|_i /2 is an approximation to {df 2 that satisfies the relation 

Aj_|_i/ 2 ((?i+i — (/j) = fi+i — fi- The term provides a second-order formulation of /j+ 1 / 2 ; various 

limiter functions may be applied to serve this function. The diffusion term in Eq. 2 is typically 
defined as follows. 

_ {uj+i + Uj) {qt+i - qj) 

2 {xi+i - Xi) 

Note that the explicit formulation of the convective flux (Eq. 3) and diffusive flux (Eq. 4) are 
independent of each other and of the source term s. The convective flux formulation exploits the 
continuity of the inviscid flux / even across shocks and slip surfaces where the dependent variable 
q may not be continuous. 

A unified algorithm for simultaneous reconstruction of the convective and diffusive flux, includ- 
ing the influence of a source term, is derived from a semi-analytic solution to the steady state form 
of Eq. 1. In this approach, the continuity of q is required in the formulation, even across shocks and 
slip surfaces which acquire finite thickness with the inclusion of physical dissipation. Semi-analytic 
reconstruction is based on a solution to Eq. 5 


dq 

dx 


d ( dq\ 

5^ +* 


( 5 ) 


where A = ^. 
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A. Constant Wave Speed Approximation (CWS) 

1. Nominal Algorithm 


Approximating the source term s, difFusivity u, and wave speed A in Eq. 5 as constants across the 
interval and introducing, for notational convenience, xl = Xi, xq = Xj_|_i/ 2 ) and xr = Xj+i 

yields the following equation. 


dq 
Ao^ - 


d?q 


dx 


= + '®0 
dx^ 


( 6 ) 


The analytic solution to Eq. 6 with constant coefficients can be derived by inspection or through 
reference to existing textbooks. However, a coordinate transformation is introduced here for two 
reasons. Eirst, a non-constant coefficient case is more easily treated in Appendix A. Second, the 
role of integration constants qo and {dq/dx)^ are more clearly defined as are their limiting forms 
as wave speed goes to zero. Consequently, the analytic solution to Eq. 6 is derived by introducing 
the following map from x to r]. 


r?(x) 

rx—xo 

= Jo exp 

Ao(a;') 


= g(«p 

Ao( 

ai-a:o) 

^0 

drj 

= exp 

1 Ao(x-a;o) 

dx 


^0 

d?r) 


Ao d-V 

dx‘^ 


I'd 

dx 


dx' 
- 1 


( 7 ) 


In the limit of Aq ^ 0 the mapping defined by Eq. 7 becomes 


rj = X — xo 

^ — 1 


dx 
d?r} ^ 
dx“^ 


( 8 ) 


0 


Substitution of Eq. 7 (or Eq. 8 if Aq = 0) into Eq. 6 yields the following. 

2 


The general solution to Eq. 9 is 

Qgeneral 


d?q 

dr]'^ 


= qo + 


so 

t'o 


dx Y 

dr]) 


(9) 


S)„" 

90+ (t)o^ 

90+ (' 


0 

exp 


(10) 


\ o { x - xq ) 

I'D 


- 1 


The particular solution to Eq. 9 is derived as follows. Eirst, integrate both sides of Eq. 9, returning 
to independent variable x on the left hand side to simplify the integration. 


J particular 


_ .sq rv / dx 
1^0 Jo \dri 

_ _so rx—xo 

1^0 Jo 


drj' 


exp 


2Ao(xQ 

I'D 


— £0 
Ac 


exp 


Ao(x-a;o) 

I'D 




( 11 ) 


- 1 


Note that f 

/ particular 


= 0; consequently, the gradient is fully defined by the general 

solution for this case where s, A , and n are treated as a constants over the interval. Integrate both 
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sides of Eq. 11 in the same manner to complete the derivation of the particular solution. 

Qparticular 


= rv 

0 \dr] J 


— 10 r ^-^0 

An J 0 


A() JO 
Aq J 0 


exp 


dr]' 


— io r^-^0 
An JC 


1 — exp 


_ 10 
An 


I'O 

Ap(a:') 


- 1 1 


VO 
An 

= ®o) - v{x)] 


{x - Xo) - ^ 


dx' 

Ao(x-xq) 

I'D 


( 12 ) 


- 1 


In the limit as Aq 


0 the particular solution in Eq. 12 goes to the expected result 

— ^0 ^ _ n2 

Qparticular — „ \X Xq) . 

ZVq 


The solution to Eq. 6 in the interval [xl,xr] may now be expressed. 

q{x Xq) Qgeneral T Qparticular 

= 10+ {'^)^d{x) + ^[{x-xo) -v{x)] 


(13) 


(14) 


The values of qq and needed for reconstruction can be solved implicitly by evaluating Eq. 14 

at the endpoints of the interval [xl, xr] using known values qi and qr. 

dq 
dx 
dq 
dx 


QR = Qo + 


m + 


Qo = 


VR - VL 


dq 

dx 


Aq 


VL - VR 


^ [{xL - xo) - Vl] 

(15) 

^ [{XR - Xo) - Vr] 

Ao 

(16) 

1 

O p 
1 

o 

(17) 

(QR- ^(xr-Vr)) 

(18) 


Note that Eqs. 17 and 18 for qo and {dq/dx)^ respectively are derived from Eqs. 15 and 16. Equa- 
tion 17 is an implicit relation in that rjR and t]r (and possibly sq) are functions of qo through Aq 
and 1 ^ 0 - A Newton iteration to solve Eq. 17 for qo is constructed as follows. 


Q = QO 


VR i^L-Xo)'^ -VL {<1R-^ {XR-Xo)'^ 


dq 

= 1 - 

AL / 

VR ^ 

OR , ' 

dqo 

dXo 1 

VR+VL j 

1 dqo P 



AL ( 

Vl ' 

1 dAo / 


dXo ( 

, Vr+Vl ) 

1 dqo p 



VR 


^o)i 


Vr-Vl ^ X 



VL 





Vr-Vl ^ Xi- 


Ao' 


(19) 




= Qo 


g 

dq 

dqQ 


The solution of Eq. 18 can be deferred until a converged value for qo is obtained from Eq. 19. The 
reconstructed interface can now be computed. 


(/ - h)o = /(go) 


( 20 ) 
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2. Discretization Considerations 


The transformed variable r] from Eq. 7 is a function of Reynolds number Rex = • 

{x - Xo) 


V = 


Rex 


(exp [Rex] — 1) 


(21) 


The cell Reynolds numbers {Rexi and Rexpt) provide a metric for resolution of diffusion dominated 
structures (boundary layers, shear layers, viscous shocks) across or within a cell. In most circum- 
stances, the cell Reynolds number is of order 1 to 10 across boundary layers and shear layers but 
can grow many orders of magnitude larger in “inviscid flow” domains (ie. regions where viscosity 
effects on the flow structure are vanishingly small). When Rexi ^ 1 (wave speed negative) and 
RexR <C — 1 the value of go from Eq. 17 is defined by the upstream conditions at xr. Conversely, 
when RexR ^ 1 (wave speed positive) and RexR <C — 1 the value of go from Eq. 17 is defined by 
the upstream conditions at xr. In numerical implementations the maximum value of \Rex\ is set 
to 100 to eliminate occurrences of numeric overflow or underflow in the evaluation of r/. 

A special case must be considered when the converged wave speed Ao = 0, even beyond im- 
plementation of the proper limiting forms in Eqs. 8 and 13. A shock is formed when > 0 and 
Xr < 0. If the shock is not resolved {\RexR\ > 1 and |i?exj^| > 1) and if Aq = 0 then the interface 


value of go is a linear average of the left and right states and the gradient is a divided dif- 

ference between the right and left states. This resulting approximation is extremely poor because 

I which is 


the transition from qr to qi through go takes place over a distance \xrs — xls\ ^ 
a factor \Re{xR — xl)\ less than the distance \xr — xr\. Without any special consideration of this 
shortcoming, the limiting form of the expressions for go (Eq. 17) and (Eq. 18) becomes 


qo = 


{XR - Xo) {qi - j^{xL - xo)^) - (xL - Xo) {qR - j^{xR - xo)^) 


XR - XL 


(22) 



0 


{(IL - ^(XL - Xo)2) - (^qR - ^{XR - Xo)2) 
XL - XR 


(23) 


The expression for in Eq. 23 may be modified by replacing values of xl and xr with xls 

and xrs respectively to account for the thin transition domain. 


^Rs ^Ls 

_ t'o 

\Xr — Xl\ 

(24) 


(xRs - XLs) 

(25) 

XLs 

II 

H 

o 

1 

to 


{xLs - XLs) 

(26) 

XRs 

= X0 + 


dq 

dx 


{(IL - ^{XLS - xo)^) - (g/? - ^{XRS - Xo)‘ 


^Ls ^Rs 


(27) 


Note that this special case arises because the wave speed at xo converged exactly to 0. Had Ao been 
initialized with a non-zero value resulting in a large value for either rjR or r]L then the converged 
interface value would approach the left or right state, respectively (see Eq. 17). In either case, 
the high gradient transition zone would occur at the opposite end of the interval, away from xo. 
This observation suggests related formulations for the analytic solution that avoids this difficulty. 
Some of these, dealing with attempts to better resolve wave speed across the cell, are presented in 
Appendix A. The next formulation is a simpler approach. 
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B. End-Point Wave Speed Approximation (EPWS) 

1. Nominal Algorithm 


In the constant wave speed (CWS) approximation of Sec. A the analytic solution to a single lin- 
earized equation centered at interface xq (Eq. 6) forms the basis for reconstruction of the inviscid 
and viscous flux. The reconstruction algorithm required a Newton iteration to accommodate the 
implicit dependence of wave speed Aq on q^. It also required a modification to account for a special 
case in which the wave speed equals zero within a shock. Both of these difficulties can be overcome 
if we reformulate on the basis of two linearized equations referenced to the right and left states 
of the interval. The origin is retained at the interface in this End-Point Wave Speed (EPWS) 
formulation. 


^4 - <-> 


The source terms carry a dual subscript LO and i?0 to acknowledge that algebraic contributions 
are taken from the appropriate right or left state while gradient contributions will carry a smaller 
stencil if evaluated at the interface xq. Transformed coordinates a{x) and j3{x) are introduced to 
take on the same role as r/(x) used in the constant wave speed approximation (Eq. 7). 


a{x) 


fX-Xo 

Jo 


exp 


I'L 



\l{x-xo) 



P{x) 


fX-xo 

Jo 


exp 


^r(x') 

‘'R 


RR 

^R 



Xr(x-Xq) 

rr 


da 


\l(x-xq) 



\r(x-Xo) 

dx 

— exp 


dx 

— exp 

fR 


d?a \l da 

dx‘^ vl dx 


d? (3 da 

dx'^ Ur dx 



(30) 


Evaluation of limiting forms as Xl ^ 0 or Ar ^ 0 as well as the general and particular solutions 
to Eqs. 28 and 29 proceed exactly as in the previous section. The left and right state reference 
solutions are: 


Qo,a — 
qo,p = 



(XR 

{qL - "^{XL-Xo)^^ 

I - 

[qR~ "^{XR- Xo)) 

(Jr 

Ot} 

(qL - ^(a^L - xo)] 

O' L 

1 — (Jl 

{qR- ^(XR-Xo)) 

( 

(Jl 

qL - ^{XL- «l)] 

i - (Jl 

- ( 

qR - ^{XR - (Xr)) 

( 

Oil 

qL - ^{XL- (Jl)) 

i — OIR 

- ( 

qR - ^{XR- Pr)) 


I3l - (Jr 


If Ai = 0 then: 


Qo,a 



0,a 


{XR - Xo) {qi - ^{XL - a:o)^) - {xr - xq) (^qr - ^{xr - xq)^) 

XR - XL 

{qL - I^{XL - Xof^ - {qR- ^{XR-Xo)^^ 

XL - XR 


(31) 

(32) 

(33) 

(34) 

(35) 

(36) 
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If Ar = 0 then: 


90,/3 


dq 

dx 


0,/3 


{XR - Xo) (^QL - ^(XL - Xo)^^ - (xL - Xo) (^QR ~ ^(xr ~ Xq)^^ 

XR - XL 

(<^L - §i^(xL - Xof^ - [qR- ^{XR- Xof^ 

XL - XR 


(37) 

(38) 


In the usual case with xq = {xl + xr)/2 expressions for qq and 
First define local cell Reynolds numbers. 



are simplified as follows. 


Rcl 

Rcr 


Rl 

\r(xl-xr) 

^R 


(39) 


Next, define reconstruction factors as a function of local cell Reynolds numbers. 


rjMO 

TfMl 


1 



exp {ReM) — i 


(40) 


with subscript M referring to the left or right node as needed. The interface values as approximated 
from the left and right states are now more simply expressed. 


qo,a = QL + {QR - Ql) VfLO + (s “ ’’/m) 

qo,/3 = qR + {qL - qR) rjRo + _ rfRo) 


(41) 


The reconstruction factors are bounded between 0 and l.The limiting forms for qq and 
as A ^ 0 for either node are written 


dx J j 


qo,a — 

\{qR + 

qL) - 

slo(xr-xl)^ 

Sul 

for 

Al - 

0 

qo,i3 = 

\{qR + 

qL) - 

SRoixR-XLr 

Sun 

for 

Ar - 

0 

d<l\ — 

1 1 



for 

Ar - 

0 

dQ\ — 

'^^Jo,l3 

qR-QL 

xr-xl 



for 

Ar - 

0 


(42) 


An arithmetic average of the left and right reference state solution for flux is defined. 


/o 

ho 


f{qo,a) + /(go,d) 
2 



+ 


2 



(43) 

(44) 


Equations 43 and 44 explicitly define the interface flux and do not require the Newton iterations 
required in the constant wave speed formulation of the previous section. 
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2. Discretization Considerations 


Numerical experiments described in the next section have been the primary basis for assessing nu- 
merical stability and accuracy. In these experiments, all implemented with implicit line relaxation, 
the simplest method for insuring stability is to guarantee that the local cell Reynolds numbers 
never exceed a specified maximum value. This limit is implemented by first defining a reference 
cell Reynolds number for the interval. 

^ max{\XL\,\XR\){\xR - xl\) 

mm{uL,nR) 

A maximum value for the reference cell Reynolds number Rcmax is specified and diffusivity at the 
left and right states is reset. 




i^R 


vr max(l, 
r'ijmax(l. 


(j)Re 

R(^max 

4>Re 

R^m.ax 


(46) 

(47) 


where 0 is a limiter bounded by 0 (no artificial diffusivity) and 1 (maximum artificial diffusivity). 
Numerical tests show that Rcmax = 8 does not engage the limiter for any reasonably resolved 
viscous layer and will result in a near perfect jump across the cell if the layer is unresolved. The 
limiting enhances stability for cases of poorly resolved boundary layers and shocks. 

Artificially increasing the diffusivity to constrain the reference cell Reynolds number from ex- 
ceeding Rcmax everywhere {(f) = 1) provides acceptable solutions to Burger’s equation. However, 
such an algorithm will undoubtedly corrupt a Navier-Stokes simulation with moderate to high vehi- 
cle Reynolds number because there is no possibility to maintain a cell Reynolds number less than 8 
across the entire boundary layer without introducing an unacceptably large number of mesh points. 
The net effect would act much the same as adding turbulent diffusivity to the solution, thickening 
the boundary layer and moderating the gradients. Numerical experiments using line-implicit relax- 
ation (next section) have been conducted to test algorithms designed to circumvent this problem. 
The tests reveal that restrictions on the cell Reynolds number are only required in the vicinity of 
a sign change in the eigenvalue (wave speed) from the left to right node. In the case of diverging 
waves (both directed away from the interface) a cell Reynolds number less than 8 is required to 
prevent the formation of an expansion shock - dissipation is the only process that connects the 
left and right states. In the case of converging waves (both directed toward the interface) the cell 
Reynolds number restriction is required to maintain stability and prevent formation of overshoots 
and undershoots across a shock. 

Conceptually, one would like to engage the requisite artificial viscosity only at interfaces where 
the product XrXr < 0. Numerical tests show that a ramp in artificial diffusivity factor is required 
from 0 to 1 in the neighborhood of the change in eigenvalue in order to achieve convergence. The 
neighborhood of a change in eigenvalue is defined as any interval in which the product (A^ -|-ti;(Ai;, — 
^r)){^r + (jj{Xr — Xr)) < 0. The factors represent a linearly extrapolated eigenvalue located to cells 
to the left and right of the respective end points. A value a; = 1 is recommended. Larger values of 
(jj tend to diffuse the shock over more cells but have beneficial effects on stability. The limiter cj) 
for engaging extra diffusivity in the neighborhood of a change in eigenvalue is now defined 


(/> = 3(/>^ - 

(p = max 


2(/>3 

0, (((j2 + uj)-{l + 2uj + 2u;2) 


ArTf 


/ (2(u^ -|- 2(jj -|- ^) 


(48) 


The variable 4> in Eq. 48 defines a normalized parameter varying between 0 and 1. A value of 
0 indicates the interval is not in the neighborhood of a change in eigenvalue or that one of the 
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extrapolated eigenvalues equals zero. A value of 1 indicates that the eigenvalues at the endpoints 
are equal magnitude but of opposite sign. The variable 4> in Eq. 48 is defined such that (p{0) = 0, 
m = 1, g(0) = 0, and *{1) = 0. 

Large values of |i?eM| cause machine zero values of reconstruction factor VfMi in Eq. 40 which 
in turn cause machine zero values of in Eq. 41 in the absence of source terms. Stability of 

the line relaxation is enhanced in these cases by resetting VfMi to 1 in Eq. 41 in order to recover 
a conventional central difference representation of diffusive gradients across neighborhood of the 
change in eigenvalue. The transition from the semi-analytic gradient reconstruction to conventional 
central difference is implemented with the following change to the reconstruction factor rjMi- 

rfMl,blended = rfMli^ — 4>) + 4> (49) 

This blended reconstruction factor may be used in place of rfMi to evaluate the gradients in Eq. 41. 


C. Application - Burger’s Equation 

Burger’s equation, the simplest test including effects of both a diffusion term and a non-linear 
convective flux term, is defined in Eq. 50. 


du 1 du^ d'^u 

dt 2 dx ^ dx"^ 


(50) 


In terms of the generic scalar conservation law presented in Eq. 1 note that q = u, f = ^,h = u^ 
and s = 0. Two sets of boundary conditions are tested: the first with a centered shock and the 
second with end wall boundary layers. Wave speed is defined hy X = ^ = u. 


1. Center Shoek 

The solution domain is defined across [—1,1]. Boundary conditions are fixed with n(— 1) = 1, 
u{l) = —1. The initial condition is defined u{x) = —x. Wave speeds are negative for x > 0 and 
positive for x < 0; consequently, waves carry the boundary values of u to the center of the domain 
where dissipation effects smooth the transition from the left to right states. The exact solution for 
the steady state limit is defined in Eq. 51. 

, , , , ax, , 

u(x) = atanh(— — ) (51) 

The coefficient a is a function of n defined with an implicit relation to satisfy the boundary condi- 
tions. 

atanh(-^) = 1 (52) 

2u 

Line relaxation is used with ^ = 0.1 for the numerical solution of Eq. 51 with waves converging 
to the center of the domain. 


2. End Wall Boundary Layers 

The solution domain is defined across [—1,1]. Boundary conditions are fixed with u{—l) = — 1, 
tt(l) = 1. The initial condition is defined u{x) = x. Wave speeds are positive for x > 0 and negative 
for X < 0; consequently, waves carry the center value of a = 0 to the boundaries of the domain 
where dissipation effects smooth the transition to the end states. The exact solution for the steady 
state limit is defined in Eq. 53. 

, , ax 

u{x) = atan( — ) (53) 
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The coefficient a is a function of v defined with an implicit relation to satisfy the boundary condi- 
tions. 

atan(^) = 1 (54) 

A constraint, 0 < a < vn/, must be imposed to avoid solutions with discontinuities. Line relaxation 
is used with ^ = 0.1 for the first 10 global relaxation steps of the numerical solution of Eq. 51 with 
waves diverging from the center of the domain. After iteration 10 the steady form of the equation 
is relaxed = 0). The converged solution takes a long time to set up as — > 0 because wave 

speed goes to zero over the domain between the end wall boundary layers. 

3. Comparisons - Shock 

The numerical solution to Eq. 50 for the center shock boundary conditions using the constant wave 
speed (CWS) reconstruction is presented as a solid black line in Eig. 1(a) for u = 1. and in Eig. 1(b) 
for ly = 0.01. There is no discernible difference between the results of the CWS formulation and 
the end point wave speed (EPWS) for zz = 1 on the scale of Eig. 1(a). The EPWS reconstruction 
uses (j) = 1 (cell Reynolds number limiting engaged everywhere) unless otherwise noted. The shock 
front sharpens as z^ — > 0. 

The broken lines show the difference between the semi-analytic formulations u and the exact 
solution Ue (Eq. 51) and between the reference formulation using Roes scheme with Yee’s Symmetric 
Total Variation Diminishing (STVD) URoe and Ue- Eor v > 0.01 all formulations are in good 
agreement with the exact solution but the semi-analytic formulations are closer to the exact solution. 
(This trend will be confirmed in subsequent tables.) 




X 


(a) u = 1.0 (b) V = 0.01, note x scale change 

Figure 1. Solution and error distribution for Burgers equation with 129 cells. 

A larger view of the solution quality for a wide range of v and cell number are summarized in 
Tables 1-6. The error norm in column 3 is defined as 

cells 

\u — Ue\ = \ui — Ue,i|Ax (55) 

i=2 

The error norm in column 4 is defined similarly for the Roe formulation. Note that the Roe scheme 
error norm in column 4 is constant for the odd number of cells with small value of z/q and cell 
Reynolds numbers greater than 100 in column 8. This result was double-checked and is caused by 
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the rounding of the shock front at the single point to the right and left of the zero point such that 
the integrated norm is a constant. Column 5 shows the number of iterations required to converge 
the residual norm of the semi-analytic formulation to the value shown in column 6. Column 7 shows 
the residual norm achieved by the reference Roe formulation in the same number of iterations. The 
last column shows the maximum cell Reynolds number achieved in any cell to provide a metric for 
the quality of resolution. 

Tables 1-3 contain error norm results and convergence data for the CWS and EPWS for- 
mulations for an odd number of cells spanning the domain. Tables 4-6 contain the same data 
for an even number of cells. Dependent variables Ui are situated at the end points of each cell; 
consequently, the odd number cell tests have adjacent nodes that span the shock for u ^ 0. The 
even cell cases have a node in the center of the discontinuity. Both instances are recorded here 
to test effectiveness of the choice Rcmax = 8 used in Eq. 46 under extremes in shock jump end 
states. The CWS and EPWS formulations show significant improvement over the reference Roe 
formulation for all simulations with v > 0.01. The Roe’s scheme has a cleaner shock jump than the 
CWS formulation for the odd number of cells and the EPWS formulation for even number of cells. 

Note that Tables 3 and 6 contain the results for EPWS where the limiter on cell Reynolds 
number is only engaged in the vicinity of a change in sign of the eigenvalue across an interval. In 
this case tc = 1 and Remax = 4. The value of 4 comes from numerical experiments indicating a 
good robustness under line relaxation without introduction of excessive smearing of shocks. The 
error norms between the two versions of EPWS are equivalent when the cell Reynolds number 
(right column) is less than 4. The globally engaged limiter (0 = 1) gives a better result than the 
restricted limiter; however, the globally engaged limiter employed a value of Remax = 8 that gives 
a near perfect shock jump for the odd number cell cases. In all other tests involving shocks with 
even number of cells or boundary layers (next section) the restricted limiter gives better results 
than the global limiter in any case where the cell Reynolds number exceeds 4. 

4- Comparisons - Boundary Layer 

The numerical solution to Eq. 50 for the end wall boundary layers using the constant wave speed 
(CWS) reconstruction is presented as a solid black line in Eig. 2(a) for i/ = 1.0 and in Eig. 2(b) for 
n = 0.01. There is no discernible difference between the results of the CWS formulation and the end 
point wave speed (EPWS) for v = 1 oi\ the scale of Pig. 2(a). Again, the EPWS reconstruction uses 
(f) = 1 (cell Reynolds number limiting engaged everywhere) unless otherwise noted. The boundary 
layer steepens as z/ ^ 0. 

The broken lines show the difference between the semi-analytic formulations u and the exact 
solution Ue (Eq. 53) and between the reference formulation using Roe’s scheme with Yee’s Symmetric 
Total Variation Diminishing (STVD) URoe and Ue- The maximum error for the v = 1.0 case 
occurs away from the wall when the boundary layer is well resolved. In contrast, the maximum 
error occurs at the wall for v = 0.01 as resolution of the high gradient region is decreased. Por 
u > 0.01 all formulations are in good agreement with the exact solution but the semi-analytic 
formulations are closer to the exact solution. (This trend is confirmed in subsequent tables.) In 
Pig. 2(b) with u = 0.01 the semi-analytic formulations are still closer to the exact solution than 
the reference formulation. The exact solution for the gradient at the half cell location next to the 
boundary re = —1 in Pig. 2(b) is (^)g = 25.996. The semi-analytic values for the gradient are 
{^)cws ~ 26.372 (1.44% error) and = 26.069 (0.28% error). The central difference 

approximation used in the reference formulation gives = 28.078 (8.0% error) . 

In contrast to the shock case there is no significant difference between odd and even number of 
cells for error norms. Tables 7 for CWS reconstruction and 8 for the EPWS reconstruction show 
similar quality of agreement with the exact solution as observed in the shock solutions. Only for 
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cell Reynolds numbers on the order of 10 to 1,000 does the reference scheme outperform the EPWS 
formulation. At Reynolds numbers higher then 2,000 the convergence of the reference scheme is 
unpredictable because no eigenvalue limiting (entropy fix) is applied. The Roe scheme diverged for 
some of the highest reynolds numbers because no eigenvalue limiting was engaged across the sonic 
line and the natural viscosity was not sufficient to prevent formation of a non-physical rarefaction 
shock. 

Table 9 contains the results for EPWS where the limiter on cell Reynolds number is only 
engaged in the vicinity of a change in sign of the eigenvalue across an interval. In this case u) = 1 
and Remax = 4 as discussed in the previous section. The error norms between the two versions of 
EPWS are equivalent when the cell Reynolds number (right column) is less than 4. The restricted 
limiter gives better results than the global limiter in any case where the cell Reynolds number 
exceeds 4. 




(a) V — 1.0 


(b) V = 0.01 


Figure 2. Solution and error distribution for Burgers equation (boundary layer) with 129 cells. 
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IV. One-Dimensional Vector Eqnation 


A one-dimensional, vector conservation law including both convection and diffusion is expressed 
by the following equation. 

do d 

A numerical solution to Eq. 56 written in conservative form must guarantee that the convective 
and diffusive flux exiting cell i is exactly balanced by that entering cell i -|- 1. Such a balance is 
enforced in the following difference equation. 


.n+l 


At 


-f 


(f - - (f - h),_ 


1/2 


(a^i+1/2 “ Xi_i/2) 


= S; 


(57) 


An analytic solution to the linearized form of Eq. 56 will be used to define the convective and 
diffusive flux at cell interfaces i ± ^ in Eq. 57. 


(9q 

di 


. dq ~ 

+ Ao^ = vqdq 
ox 


dx'^ 


+ So 


(58) 


where q is the conserved variable set and q is the set of primitive variables defining the diffusion 
flux where Tdq = dq. Different strategies have been explored to diagonalize this system with 
various choices of artificial diffusion to remove the singularity from B and based on various choices 
for the primitive variable set q. A characteristic variable formulation provides a foundation that 
maintains the physical wave speed and diagonal dominance (diagonal elements are all positive 
definite, off-diagonal elements are negative) in the linearization of the diffusion flux. 


(9q 

~di 


5q clq 

dt ^ dx 
dq 




z/oBT“^T 

1^0 (bT-1 


d^q 

dx"^ 


+ So 


d^q 
0 dx'^ 


+ So 


5q' 


, * dq' 

+ Ao^ = t'o 
ox 


R_igT-iR 


R-l 

^ 9 -T -Iv- OSQ 

0 


(59) 


The steady form of Eq. 59 is approximated as follows with T o equal to the diagonal of 
to simplify derivation of an analytic solution. 


(^r-^bt-ir) 


0 


Ao^ =i/oTo^+R-Vo + ^^o(R“'BT-1R-T 
ax ax^ 


d^\ 
dx^ y , 


(60) 


The last term of Eq. 60 may be combined with the source term. 


A. End-Point Wave Speed (EPWS) - Nominal Algorithm for Systems 

The EPWS formulation of a semi- analytic solution will be developed here for systems of equations 
because it yields an explicitly defined interface. The solution expanded with respect to the left end 
point condition is developed in detail. The solution with respect to the right end point condition 
may then be inferred. The nomenclature closely follows that used in Sec. Ill - B. 

The analytic solution to Eq. 60 using left state reference conditions is expressed. 

(61) 


16 of 52 


American Institute of Aeronautics and Astronautics 



The elements of diagonal matrix T'q are defined as a function of the jth element of A/,, Xlj, and 
the jth element of r-iX l, 


'4^a,j — 




^L,j{x-Xo) 

''L,j 
X — Xq 


- 1 j Xi,j / 0 

Alj = 0 

The elements of diagonal matrix are defined 

X _ / 1Vl7 ~ ~ ^L,j 7^ 0 

’'“I >^L,j=0 

Conversion from characteristic variables to conserved variables follows. 


(62) 


(63) 


RZ'(q«-qa,o) = +^a^l^S0L (64) 

qo = + Rl^qRj^^sql (65) 

The left and right values of qa are functions of T'o(a^L) = ^^a,L, '^aixR) = ^a(xL) = ^a,L 

and ^a{xR) = ^a,R- 


qL = qa,0 + + Rl^o,lRl ^SqL 

q/? = qa,0 + RL^a,i?R hS).* RL^a,ijR2.^®0L 

The interface conditions for q^^o and can be solved from Eqs. 66 and 67. 


( ^ 

I dx 


a,0 


Rl (^a,L - ^a,R) ^ 

[R-Z^ (qi- “ q/j) “ (^a,L - ^o,r) RZ^sql] 


( 66 ) 

(67) 


( 68 ) 


q «,0 = Rl (^a,i? — ^o,l) ^ 

[^o.flRZ^qi “ ^a.iRZ^q^? “ ('^a,R^a,L ~ ^a,L^a,i?) RZ^^Ol] 


(69) 


The right end point expansion follows the same pattern as the left. The elements of diagonal matrix 
^/3 are defined 




^7exp[AtM0^ 



X — Xq 


Xrj / 0 
Xrj = 0 


(70) 


The elements of diagonal matrix are defined 


^13, j 


[{x - xo) - ^ 0 ,j] Xrj / 0 

-27ii-(^-^o)^ Arj=0 


The interface conditions for q^ q and 



are defined. 


dx J 


13,0 


Rfl (^/3,L - ^/3,i?) ^ 

[RZ^ (q^ “ q«) “ i^f3,L - ^(3 ,r) RZ^SOij] 


(71) 


(72) 
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( 73 ) 


Q/3,0 — R--R ^ 

- ^/3,LR-R^q_R - (^/3,R^/3,L - ^/3,L^/3 ,k) R-r^SQ/j] 

Additional simplifications expressed below (equivalent to scalar equation Eq. 41) can be intro- 
duced if one assumes xq = ^{xr + xr). 

q«,o = qL + RlD^/loRl^ (q/j — qi) + (|l - Dr/Lo) (iCij — xl) Sol 

q/3,0 = qi? + RfiDr/LoR)?^ (qL - q_R) + RrA];j^ ( 2 I - Rj?^ (a^L - a^ij) soij 

(®)„,o = RlD,;iiRj‘^+RiAj'(I-D,,„)RZ‘ 80 L 

(=)«,o = R«D-/hiRh S3f +KRA;;‘(I-D,/Hi)Ri'%R 

where DrfMO and D^/mi are diagonal matrices of reconstruction factors whose elements are 
bounded between 0 and 1 (with M a dummy variable for R or L) corresponding to scalar equation 
Eq. 40 whose elements are defined 


^rfMOJ 




( ReM,i 

exp 


exp ReM,j 

exp {Reuj) ~ 1 


RSRj = 




ReRj = 


' ^L,j 

_ >^R,j(^L-XR) 

' ~ ^R,j 


Note that vmj is the jth element of diagonal matrix um"^ M and Xmj is the jth element of diagonal 
matrix of eigenvalues Am- 

The limiting forms for the jth diagonal elements when \j ^ 0 in Eq. 74 are defined 


drfMOJ — 2 — 0 

^(i-i/Roj) = for \Lj = 0 

^(l-<i./R0,,) = f” ^Rr = » 

drfMi,j = 1 for Xmj = 0 

\b- - drfMij) = 0 for Xmj = 0 


for Xrj = 0 
for Xmj = 0 

for Xmj = 0. 


B. Discretization Considerations for Systems of Equations 

A limited value of diffusivity is required to maintain a robust capture across shocks and expansions 
where an eigenvalue passes through zero. The limiter is based on a maximum value of a reference 
cell Reynolds number Rcmax- The maximum is set to 8 if the limiter is applied globally as discussed 
earlier for Eq. 46. The maximum is set equal to 4 if the limiting is introduced only in the vicinity 
of a change in sign of eigenvalue. The reference cell Reynolds number is a scalar value that is still 
defined with Eq. 45, where |A| is taken as the maximum norm of the diagonal matrix of eigenvalues. 

^ max(|AL|,|Aj;|)(|xfl-XL|) 
mm{i/L,T^R) 
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The difFusivity is reset as originally discussed in Eq. 46. 


VL ^ jvimax(l, 

UR ^ URma,x{l,^^) 


( 46 ) 


where (p is a limiter bounded by 0 (no artificial difFusivity) and 1 (maximum artificial difFusivity) 
as originally presented in Eqs. 48. Equation 48 is slightly modified below in order to emphasize 
that one is searching for a maximum over all j elements of the matrix of eigenvalues. 


(p = 3(/)2 - 2(p^ 

(^(^2 + _ (1 + 2 w + 7 ( 2^2 + 2 (u + i ) 

Finally, the transition from the semi- analytic gradient reconstruction to conventional central 
difference is implemented with the following change to the diagonal matrix of reconstruction factors 
D rfMi as presented earlier in scalar Eq. 49. 

^rf Ml, blended — (1 4^)^rfMl T 01 (^0) 

This blended reconstruction factor may be used in place of T)rfMi to evaluate the gradients in 
Eq. 74. 


= max < 0, max 



C. Application - Quasi, One-Dimensional Nozzle 

The governing equations for quasi one-dimensional nozzle flow are expressed 

where A(x) is the area distribution of the nozzle and 

pu 

f = I p -|- pu‘^ 
puH 

0 


h = 


i du 
3 ^ dx 

I 9T 

L Vr dx A 


S = 


0 

Pdx 

0 


The conserved variable vector q and primitive variable vector q are defined 



p 


p 

q = 

pu 

, q = 

u 


. PE _ 


T 


Eq. 81 may be rewritten to move all nozzle area terms to the source. 

clq . 5q ^d^q 
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(83) 


(84) 


(85) 
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where the source term is now defined by 


d log A 

s= 

[ J 

and relations for a calorically perfect gas are assumed. 

p = (5pe 
e = CyT 

P = J-l=Cp/Cy-l 

= ryp/p 

The linearization employs matrices A and B defined below. 


dq 


0 1 0 
— 1) u{2 — P) P 

uiv?^ — H) —(du^ + H Pu + u 


( 87 ) 


( 88 ) 

(89) 

(90) 

(91) 


(92) 


dh 



uB = u 


0 0 0 
0 0 
0 ^pu _ 


(93) 


Transformations from primitive to conservative variables (Tdq = dq) and conservative to primitive 
(T~^dq = dq) are defined 



1 

0 

0 


1 

0 

0 

T = 

u 

P 

0 

T-1 ^ 

u 

p 

1 

p 

0 


E 

pu 

pCv 


1 

1 

1 

4T 

1 

1 

1 

pCv - 


Diagonalization of A = RAR ^ is implemented with 



2 

1 

1 


■ P-Pi 

Pu 

-p 


2u 

u + c 

u — c 

R-1 = 

P ^ — uc 

c — Pu 

p 


_ 

H + cu H — cu 


_P^ + UC 

—c — Pu 

p 


A = 


u 0 0 

0 u + c 0 

0 0 u — c 


(94) 


(95) 


(96) 


All matrices appearing in Eq. 60 are defined for the quasi one-dimensional nozzle problem. The 
product R“^BT“^R and its diagonal T in this equation are evaluated. 


R_igT-i 



2 -P -P 
-2 P+^ p-^ 
-2 P-^ . 


r = 


1 


2 0 0 

0 ^ 0 

0 0 p+^ 


(97) 


( 98 ) 
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Three nozzle flows are simulated with area ratios of 20, 200, and 2,000 and respective exit 
pressure to stagnation pressure ratios of 0.1, 0.01, and 0.001. The gas model (MKS units) uses 
ratio of specific heats 7 = 1.4, molecular weight = 28.8, Prandtl number Pr = 0.72, and 
Sutherland’s law for viscosity, 

/i(T) = 0.00001458205T3/V(^ + HO.333). (99) 

The baseline cases specify Tgtag = 2,000. and Pstag = 10®. The baseline specifications result in 
reference cell Reynolds numbers exceeding 8 everywhere - even with the finest grid. An off baseline 
case specifies Pstag = 10 ^ which brings reference cell reynolds numbers to order 1 or less in the 
vicinity of a captured shock in the diverging section of the nozzle. 

The area distribution is defined in Fig. 3. It has zero slope at the entrance and exit to simplify 
formulation of extrapolated boundary conditions at x = ±1. (The change in characteristic variables 
across the boundaries are set to zero when associated eigenvalues direct waves out of the nozzle.) 
The exponent (4) provides a very smooth, slow variation of area at the throat even for the largest 
area ratio tests. 



Figure 3. Area distributions used in quasi-one-dimensional nozzle tests. 


At the subsonic inflow boundary 


ifi = H. 


P^Pstag 


stag 
Pstag Pi 


R3,i{P1 — P 2 ) + R3,2{{pu)l — {pu) 2 ] + Rz^ 3 {{pE)i — {pE) 2 ) — 0. 


Exit pressures are specified to cause a standing shock in the diverging portion of the nozzle. Thus, 
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at the subsonic outflow boundary 


PN 

Ri,i{pn-i — Pn) + Ri,2{{pu)n-i — {pu)n) + Ri,3{{pE)n-i — {pE)n) 
R2,i{pn-i — Pn) + R2,2{{pu)n-i — {pu)n) + R2,z{{pE)n-i — {pE)n) 


Pexit 

0 

0 . 


Implicit line relation is used to drive the residual of the discretized governing equations to zero. 


Ati 


+ 


^i+l/2 (f - ^)i+l/2 ~ ^j-1/2 (f - h)j_ 


1/2 


(^i+1/2 “ ^i- 1 / 2 ) 


= S; 


A variable Ati at each cell is governed by 


A 

Ati 


Ui\ + Ci 

CFL*^ 


(100) 


( 101 ) 


where CFL'^ = min(1000., l/(Residual"'). Typically the global residual starts at a value of order 

1 and decreases as convergence is achieved to 10“^®. Updates to the conserved variable Aq are 
preprocessed by a safety factor Cgafe- The value of Csafe has a maximum value of 0.25 but will be 
made smaller to prevent the change of positive definite quantities p or e from exceeding 90% of 
their current value. Larger, more aggressive settings for Cgafe will work for some of the smaller area 
ratio problems but are not documented here. 

Figures 4-6 present the Mach number distribution and error in Mach number from inflow to 
shock (isentropic domain) for each nozzle. Solid lines indicate results from Roe’s scheme^ with a 
Symmetric Total Variation Diminishing (STVD)^^’^^ limiter. Dashed lines indicate semi-analytic 
reconstruction using the end-point wave scheme (EPWS). Variations of EPWS with cj) = 1 (global 
application cell Reynolds number limiting) and cj) = with cj = 4 (local application of cell 

Reynolds number limiting within 4 cells of a critical point) are tested. Results for five levels of grid 
refinement are documented. The global view of Mach number shows little difference between either 
of the EPWS formulations and Roe’s method. All of the algorithms remain stable from nearly 
incompressible to hypersonic domains. Significant differences are evident as scale is refined. The 
exact solution in Pigs. 4 - 6 refers to an inviscid system of equations. All of the methods produce 
larger errors in the throat (transonic domain) than in either the subsonic or supersonic domain 
- probably associated with the introduction of the limiting algorithm. Error magnitudes increase 
with increasing maximum area ratio as a result of larger rate of area variation. 

Error at x = —0.5 (subsonic domain) for the area ratio 2,000 case is plotted as a function of 
mesh size for the Roe and EPWS schemes in Pig. 7. The slope of the error curve is approximately 

2 for Roe and EPWS with cj) = (p{uj) indicating second-order accuracy. The slope of the error 
curve for EPWS with 0 = 1 is approximately 1 indicating first-order accuracy. These results 
are representative of those observed for other area ratios and at other positions along the nozzle 
away from the throat and the shock. The constant cj) formulation adds a dissipation term linearly 
proportional to the cell Reynolds number up to the point where the cell Reynolds number is less 
than 8. Por the present test case where the cell Reynolds number exceeds 8 everywhere the added 
diffusivity works like a first-order dissipation. The Roe’s scheme engages a limiter and the EPWS 
scheme with (p = 4>{cv) engages added diffusivity only at isolated locations in the domain so that 
second-order properties appear to recover away from these isolated locations. The EPWS scheme 
achieves this order property with a 3 point stencil while the Roe scheme requires a 5 point stencil 
for the STVD limiter. All subsequent figures for the nozzle problem will use the EPWS scheme 
with (j) = 

Pigure 8 presents temperature distributions in the nozzle for the baseline stagnation pressure 
{Pstag = 10® Pa) and a lower stagnation pressure {pstag = 10^ Pa). The cell Reynolds number 
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(a) Mach number, Roe - solid, EPWS 0 = 1- dashed (b) (M — Mf,xact)/Mf.xact, Roe - solid, EPWS 0=1- 

dashed 




(c) Mach number, Roe - solid, EPWS 0 = 0(w) - dashed (d) {M — Mexact) / M^xact, Roe - solid, EPWS 0 = 0(o;) - 

dashed 


Figure 4. Mach number distribution and error for = 20. 

^throat 


crossing the shock is of order 1000 for the baseline stagnation pressure (Fig. 9(a)) and is of order 1 
for the lower stagnation pressure (Fig. 9(b)). Shock structure is resolved for cell Reynolds numbers 
of order 1 or smaller while an abrupt jump spanning two cells is obtained for the higher cell Reynolds 
numbers (limited to a value of 8 in Eq. 45 and Eq. 46 with <0 = 1). In both cases, profile convergence 
is observed with grid refinement. Also, the Roe scheme (solid line) and the semi- analytic EPWS 
scheme (dashed line) converge with grid refinement. 

Details of the temperature distribution at the throat and behind the shock are presented in 
Eigures 10 - 12. Both the EPWS and Roe’s scheme show a temperature overshoot approaching the 
throat for the two largest area ratios with the largest source terms. The overshoot is less pronounced 
for the EPWS reconstruction which includes the source term in the reconstruction process. Both 
the Roe and the EPWS schemes present a smooth transition across the shock. 

The EPWS scheme requires approximately four times the CPU time per relaxation step for the 
coupled viscous and inviscid flux reconstruction as compared to the Roe scheme as seen in Table 10. 
Only results for the maximum area ratio equal to 20 are shown because the limiter used in the 
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(a) Mach number, Roe - solid, EPWS 0 = 1- dashed (b) (M — Mf,xact)/Mf.xact, Roe - solid, EPWS 0=1- 

dashed 




(c) Mach number, Roe - solid, EPWS 0 = 0(w) - dashed (d) {M — Mexact) / M^xact, Roe - solid, EPWS 0 = 0(o;) - 

dashed 


Figure 5. Mach number distribution and error for = 200. 

^throat 


Roe algorithm hung for the larger area ratio cases. The EPWS scheme usually converges in fewer 
relaxation steps (column 5) so that total time to convergence (column 3) is roughly equivalent - a 
result somewhat dependent on the magnitude of ringing and time to damp same associated with the 
limiter. Both the EPWS and Roe reconstruction use identical relaxation algorithm. The coarsest 
grid (36 cells) is started from a crude initialization; subsequent grids are initialized by interpolation 
from the previous converged solution. 
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(a) Mach number, Roe - solid, EPWS 0 = 1- dashed (b) (M — Mexact)/Mexact, Roe - solid, EPWS 0=1- 

dashed 




(c) Mach number, Roe - solid, EPWS 0 = 0(w) - dashed (d) {M — Mexact) / M^xact, Roe - solid, EPWS 0 = 0(o;) - 

dashed 


Figure 6. 


Mach number distribution and error for 

^throat 


2 , 000 . 
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Figure 7. Grid convergence at a; = —0.5 for the area ratio 2000 case indicating order properties of the 
EPWS and Roe formnlations. 



(a) petag ~ 10® Pa, Roe - solid, EPWS (f) = 4i(ui) 


(b) pstag = 10® Pa, Roe - solid, EPWS </> = </>(w) - 


dashed 


dashed 


Fignre 8. Temperature distribution for = 2, 000. 

■^throat 
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(a) 

pstag - — ■ 10® Pa, Roe - solid, EPWS - dashed (b) pstag = 10® Pa, Roe - solid, EPWS - dashed 


Figure 9. Cell Reynolds number distribution for = 2, 000. 

■^throat 




(a) Throat, Roe vs EPWS cf> = 


(b) Shock, Roe vs EPWS (j) = 4 >{uj) 


Figure 10. 


Temperature distribution for 

■^throat 


2 , 000 . 
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36 




(a) Throat, Roe vs EPWS cj) = cf>{uj) (b) Shock, Roe vs EPWS (j) = 4 >{uj) 


Figure 11. Temperature distribution for = 200. 

^throat 


36 




(a) Throat, Roe vs EPWS 0=1 with 0 = 0(tj) (b) Shock, Roe vs EPWS 0 = 0(w) 


Figure 12. Temperature distribution for ^ = 20. 
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V. Multi-Dimensional Vector Equation 


Equations are written in a locally convenient, Cartesian coordinate system (X,Y,Z). In the 
applications described herein X is aligned with a line segment xr — xr with origin at the midpoint 
of the line segment. In future work, it may be advantageous to align X with convective velocity 
u with origin at the centroid of a tetrahedron. In either case, the transformation from reference 
system x to local system X and back is given by: 


’ dX ' 


1 

1 


dx 

dY 

= 

lx ly Iz 


dy 

dZ 


rrix rriy rriz 


dz 


dx 


I 

H 

1 


’ dX ’ 

dy 

= 

TT-y ly 


dY 

dz 


1 

1 


dZ 


Velocities are transformed in like manner. 


’ U ’ 


'^x 


u 

V 

= 

lx ly Iz 


V 

_ w _ 


rrix rriy rriz 


w 


u 


1 

1 


’ U ’ 

V 

= 

Tiy ly ^^^y 


V 

w 


1 

1 


_ w _ 


(102) 


(103) 


(104) 


(105) 


In the case of X aligned with the line segment the transformation is defined with Ux = dx/ds, 
Uy = dy/ds, and Uz = dzjds where ds = dx"^ + dy^ + dz^. Unit vectors I and m are mutually 
orthogonal. 

The conservation equations include an inviscid part f, a viscous part h and a source term s 
that govern the evolution of conserved variable q. 


5q 5(fi-hi) a(f2-li2) ^(fs-ha) 
dt^ dX ^ dY ^ dZ 


(106) 


A solution to a steady, linearized approximation to Eq. 106 is sought of the form q(V, Y, Z) = 
qo+qi(V)+q2(V)+q3(Z). The viscous terms defining h involve derivatives of q with respect to X, 
Y, and Z (h[q, d^/dX, dq/dY, dc{/dZ]). Because of the separation of variables in the assumed form 
for q (and implicitly q) the only parts of hi to survive the linearization involve second derivatives 
of q with respect to X. In like manner, the linearization of h2 include only second derivatives with 
respect to Y and the linearization of ha include only second derivatives with respect to Z. All cross 
derivative terms are annihilated. The steady, linearized form of Eq. 106 thus becomes: 


where 


9qi ^ <9q2 ^ ^ 

Ai— — - ^ + A2 t^ - l^B22VX7y + Astt^ - Z^B33— ^ - So 


dX 


dX^ 


dY 


dY‘i 


'dX 


dZ^ 


A, 




dtj 

dq 

dhi 


(107) 


(108) 

( 109 ) 
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Isolating the scalar factor v = p keeps elements of Bjj of order 1. 

The semi-analytic solution for flux reconstruction must satisfy appropriate boundary conditions 
in addition to Eq. 107. However, boundary values are only defined at nodes. Consequently, there 
is an infinite degree of freedom in how one may construct a solution that needs to satisfy boundary 
conditions at isolated nodes. It is this freedom that accommodates a solution employing separation 
of variables. This same freedom introduces ambiguity in how one handles the effect of source term 
So on the solution of qi, q2, and qa. The approach taken here is to distribute the effect of sq 
evenly over the contributions from the three coordinate directions. Thus, solutions are derived for 
separable pieces of Eq. 107. 


Ai^ 

^ dX 

~ d^qi 1 

- '^^ 11 ^X 2 + 3 "° 

(110) 


- d‘^q2 1 

- ^^^22 + gSo 

(111) 

A3t® 

^dX 

- 52 % 1 

- >'B33 + 383 

(112) 


Solutions to Eqs. 110 - 112 were defined in the previous section in Eqs. 66 - 67 with RjT'j(Xj) and 
^i(Xi) defined in Eqs. 62 and 63. The matrices for three-dimensional flow of a perfect gas and for 
a mixture of thermally perfect gases in thermochemical nonequilibrium are defined in Appendix B. 


qi = qi,o + Ri^'i(X)Rr' (ll) +Ri$i(A)R^i^So (113) 

q2 = q2,o + R2^2(E)R2i + R2$2(T)R2 (114) 

qa = q3,o + R3^3(^)Ra' (^) + R3^3(^)Rs (115) 

The solution for q is obtained by summing Eqs. 113 - 115. 

q = q„ + RiM)i(X)Rr‘(^)^ + R 2 >I. 2 (y)Rj-‘(|a)^ + R3>I.3(Z)R3‘(§)^ 

+ 1 (Ri$i(A)Rr^ + R2$2(E)R2 ^ + R3$3 (^)Rs ') So 

Equation 116 requires evaluation of four vectors: qo, i (^) ’ (^) ' 

can be evaluated implicitly by applying known boundary values q^, Qb, Qc, and qn at nodes A, 
B, C, and D. 

qA = 

+ 

qB = 

+ 

qc = 

+ 

qD = 

+ 
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qo + Ri^i.cRi ^ (Ix)q + R2'®'2,cR2 
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f 8q \ 

y9z)o 


( 117 ) 
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Equation 117 can be simplified as shown below in Eq. 118, however, it still requires an implicit 
solution of a 15 X 15 matrix for even a perfect gas. Additional simplifications are possible to 
eliminate one of the remaining gradient vectors resulting in a 10 x 10 matrix for a perfect gas 
but ultimately an implicit solution is required, even for the EPWS algorithm. (In the EPWS 
formulation, one would substitute each of the nodal subscripts A, B, C, D in turn for the centroid 
subscript 0.) Consequently, the formulation of a fully three-dimensional reconstruction requires 
significantly more operations to handle the implicit solution than the one-dimensional systems 
which are handled explicitly. 

qA-qB = + 

+ 5 [Ri(^ 1,A — ^l,s)Ri ^ + R2(^2,A - ^2 ,b)R2 ^ + R3(^3,A “ ^3,b)R3 §0 

qs-qc = Rl(^'l,B-^'l,c)Rr'(|l)Q + R2('I'2,B-'I'2,c)R2“'(|^)^ 

+ R3(^'3,B-^3,c)R3'(||)q 

+ 5 [Ri(^ 1,B — ^l,c)Rr^ + R2(^2,B - ^2 ,c)R 2 ^ + R3(^3,B “ ^3,c)Rs §0 

qc-qn = Rl(^'l,c-^l,D)Rr'(S)g + R-2(^2,C-^2,D)R2"'(#)^ 

+ R3(^^3,C - ^3,d)R3"' (II)^ 

+ 5 [Ri(^ 1,C — ^l,D)Rr^ + R2(^2,C — ^2 ,d)R 2 ^ + Rsl^S.C “ ^3 ,d)R3 §0 

(118) 

dqA-B = Rld'I'l,A-BRr' (|r)o + R-2d^2,A-BR2“^ (|^)g 

+ ^ [Rid^l A-sRi ^ + R2d^2,A-sR2 ^ + Rs^^^S.A-sRs So 
dqB-C = KldAf,^B-C^T" {§^)^ + ^2dAf2,B-C^2^ 

+ R3d'I'3,S-cR3“'(||)° ° (119) 

+ g [Rid$l^B-cRl ^ + R2d^2,_B-cR2 ^ + RsC?^3,_B-cR3 ^0 

dqC-D = {§^)^ + ^2dAf2,C-D^2^ 

+ K,{dA>s,C-D^ 3 ^{^)^ 

+ ^ [Rid$l^C'_£)R3 ^ + R2(i^2,C-DR2 ^ + R3'^^3 ,C-dR3 ^0 

An explicit formulation can be recovered if an edge specific, one-dimensional (as opposed to 
volume specific, three-dimensional) reconstruction is accepted. In this case, align X between the 
nodes A and B. Independent variables Y and Z are identically equal to zero on this edge. Conse- 
quently, 'J' 2 (^)j "^ 3 (Z), $ 2 (^)) and ^ 3 (Z) are identically zero on this edge. Boundary conditions 
at nodes A and B are sufficient for defining the semi-analytic flux at any point on this line seg- 
ment. The reconstruction algorithm now plays out exactly as described in the previous section 
for one-dimensional systems with appropriate substitution of the matrices for three-dimensional 
gas flow as defined in Appendix B. Preliminary tests of this simplification (which in fact sacrifices 
the original intent of multidimensional reconstruction) exhibited stability problems thought to be 
associated with inadequately moderating the diffusivity in the vicinity of edges where eigenvalues 
change sign. 

Programming the implicit formulation at the centroid of a tetrahedron involves significant 
changes to the basic structure of existing codes (FUNSD^^ in this case). The loop over edges 
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to create inviscid flux disappears. Computation of the gradient field used for second-order correc- 
tions can also be omitted. Flux reconstruction occurs in a loop over centroids in which inviscid 
and viscous flux contributions are computed simultaneously. The complexity of the reconstruction 
algorithm should not be trivialized here. If the number of relaxation steps required to evaluate the 
semi-analytic solution at the centroid is greater than the ratio of total edges to total centroids then 
the work for semi-analytic reconstruction in 3D carries a penalty (as already documented in the 
quasi-lD nozzle tests). Application of a directionally dependent diffusivity across critical points 
must also be derived and tested. The ID nozzle solutions demonstrate that the approach can be 
used to capture strong discontinuities with good retention of accuracy. The flexibility to provide 
truly three-dimensional reconstruction as well as the unified formulation of viscous and inviscid 
flux provides motivation for continued development. 

VI. Concluding Remarks 

Semi-analytic reconstruction uses the analytic solution to a second-order, steady, ordinary differ- 
ential equation (ODE) to simultaneously evaluate the convective and diffusive flux at all interfaces 
of a finite volume formulation. Influence of source terms on the reconstructed flux is explicitly 
included in the formulation. The second-order ODE is itself a linearized approximation to the 
governing first- and second-order partial differential equation (PDE) conservation laws. Thus, 
semi-analytic reconstruction defines a family of formulations for finite volume interface fluxes using 
analytic solutions to approximating equations. 

Limiters are not applied in a conventional sense; rather, diffusivity is adjusted in the vicinity 
of changes in sign of eigenvalues in order to achieve a sufficiently small cell Reynolds number in 
the analytic formulation across critical points. Cell Reynolds numbers of order 1 are required to 
analytically span a captured shock without introduction of significant oscillations. This target cell 
Reynolds number is defined by numerical experiment in the quasi- ID nozzle test problems but 
the result may be derived analytically for Burgers equation. The target cell Reynolds number is 
achieved by raising the diffusivity only in the neighborhood of critical points but conceivably grid 
adaptation could be used to decrease cell size as well. 

The reconstruction stencil for the new algorithm is compact, much like a finite element approach. 
The current approach, however, requires a second-order PDE; it cannot be applied directly to the 
Euler equations. The approach is finite- volume-based in the sense that conservative flux formulation 
is the primary goal and reconstructed fluxes are CO continuous across cell walls. 

Several approaches for application of semi-analytic reconstruction for the solution of one- 
dimensional scalar equations are introduced. Results are compared with exact analytic solutions 
to Burger’s Equation as well as a conventional, upwind discretizations using Roe’s scheme. One 
approach, the end-point wave speed (EPWS) approximation, is further developed for more com- 
plex applications. One-dimensional vector equations are tested on a quasi one-dimensional nozzle 
application. The EPWS algorithm has a more compact difference stencil than Roe’s algorithm 
but reconstruction time is approximately a factor of four larger than for Roe. Though both are 
second-order accurate schemes. Roe’s method approaches a grid converged solution with fewer grid 
points. Work is ongoing to improve these metrics for semi- analytic reconstruction. 

Reconstruction of flux in the context of multi-dimensional, vector conservation laws including 
effects of thermochemical nonequilibrium in the Navier-Stokes equations is developed as an em- 
bedded implicit solution. Programming the implicit formulation at the centroid of a tetrahedron 
involves significant changes to the basic structure of existing codes and has not yet been completed 
for 3D cases. The ID nozzle solutions demonstrate that the approach can be used to capture strong 
discontinuities with good retention of accuracy. The flexibility to provide truly three-dimensional 
reconstruction as well as the unified formulation of viscous and inviscid flux provides motivation 
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for continued development. 

The approach taken throughout has been to linearize with respect to characteristic variable qb 
Yet review of a rich tradition of analytic solutions in boundary-layer theory (e.g Blasius solution) 
shows that transformations of both independent variables (Howarth-Dorodnitzyn) and dependent 
variables have been introduced to generate solutions. One may speculate that lessons learned in 
that tradition may be exploited to provide better reconstruction in the algorithms developed here. 
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Table 1. Constant Wave Speed - Shock - Odd Number of Cells 
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Table 2. End Point Wave Speed - Shock - Odd Number of Cells - (j> = 1 
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3.8E-10 

1.2E-01 

33 

2.8E-05 

1.4E-04 

13 

2.9E-11 

7. IE- 11 

6.1E-02 

65 

7.2E-06 

3.2E-05 

18 

2.6E-11 

3.5E-11 

3.1E-02 

129 

1.8E-06 

7.3E-06 

26 

5.4E-11 

5.9E-11 

1.6E-02 

5 

l.lE-03 

4.1E-02 

12 

2.2E-11 

l.OE-09 

4.0E-h01 

9 

5.7E-04 

3.9E-02 

13 

1.8E-11 

1.8E-05 

2.2E-h01 

17 

l.OE-03 

3.7E-02 

14 

4.4E-11 

6.8E-05 

1.2E-h01 

33 

l.OE-02 

2.5E-02 

18 

2.4E-11 

1.7E-05 

6.1E-H00 

65 

9.7E-03 

l.lE-02 

24 

2.6E-11 

6.8E-07 

3.1E-H00 

129 

1.5E-03 

3.2E-03 

32 

6.8E-11 

3.7E-09 

1.6E-h00 

5 

l.lE-03 

4.0E-04 

12 

2.2E-11 

1.4E-10 

4.0E-H03 

9 

5.8E-04 

4.0E-04 

13 

1.8E-11 

2.8E-07 

2.2E-H03 

17 

3.1E-04 

4.0E-04 

14 

4.6E-11 

4.3E-07 

1.2E-h03 

33 

1.6E-04 

4.0E-04 

18 

2.4E-11 

1.3E-06 

6.1E-H02 

65 

8.1E-05 

4.0E-04 

23 

2.9E-11 

1.3E-06 

3.1E-H02 

129 

4.1E-05 

4.0E-04 

31 

3.0E-11 

1.3E-06 

1.6E-H02 

5 

l.lE-03 

4.0E-06 

12 

2.2E-11 

3.1E-10 

4.0E-H05 

9 

5.8E-04 

4.0E-06 

13 

1.8E-11 

4.4E-10 

2.2E-H05 

17 

3.1E-04 

4.0E-06 

14 

4.6E-11 

6.6E-09 

1.2E-H05 

33 

1.6E-04 

4.0E-06 

18 

2.4E-11 

5.0E-10 

6.1E-H04 

65 

8.1E-05 

4.0E-06 

23 

2.9E-11 

5.0E-10 

3.1E-H04 

129 

4.1E-05 

4.0E-06 

31 

3.0E-11 

2.7E-09 

1.6E-h04 

5 

l.lE-03 

3.9E-08 

12 

2.2E-11 

2.4E-08 

4.0E-H07 

9 

5.8E-04 

4.0E-08 

13 

1.8E-11 

1.8E-10 

2.2E-H07 

17 

3.1E-04 

4.1E-08 

14 

4.6E-11 

7.0E-09 

1.2E-H07 

33 

1.6E-04 

4.0E-08 

18 

2.4E-11 

9.9E-10 

6.1E-H06 

65 

8.1E-05 

4.0E-08 

23 

2.9E-11 

1.3E-10 

3.1E-H06 

129 

4.1E-05 

4.0E-08 

31 

3.0E-11 

2.7E-12 

1.6E-h06 
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Table 3. 

V 

End Point Wave Speed - Shock - 

cells \u — Ue\ \uRoe — Ue\ Iter. 

Odd Number of Cells - 4> = 

L 2 L2,Roe Ax/v 


5 

9.7E-08 

9.4E-07 

3 

8.5E-11 

2.5E-10 

4.0E-03 


9 

3.3E-08 

2.8E-07 

4 

1.4E-13 

l.OE-12 

2.2E-03 

100 

17 

9.5E-09 

6.5E-08 

4 

1.2E-12 

3.0E-12 

1.2E-03 

33 

2.5E-09 

1.4E-08 

4 

1.3E-11 

1.6E-11 

6.1E-04 


65 

6.6E-10 

3.2E-09 

5 

2.3E-13 

4.7E-13 

3.1E-04 


129 

1.7E-10 

7.4E-10 

5 

8.9E-12 

9.4E-12 

1.6E-04 


5 

1.2E-03 

9.2E-03 

7 

7.9E-11 

8.7E-10 

4.0E-01 


9 

3.7E-04 

2.8E-03 

8 

7.7E-11 

4.5E-09 

2.2E-01 


17 

l.lE-04 

6.6E-04 

10 

3.7E-11 

3.8E-10 

1.2E-01 

i.OO 


33 

2.8E-05 

1.4E-04 

13 

2.9E-11 

7. IE- 11 

6.1E-02 


65 

7.2E-06 

3.2E-05 

18 

2.6E-11 

3.5E-11 

3.1E-02 


129 

1.8E-06 

7.3E-06 

26 

5.4E-11 

5.9E-11 

1.6E-02 


5 

1.6E-01 

4.1E-02 

21 

4.7E-11 

l.lE-14 

4.0E-f01 


9 

9.1E-02 

3.9E-02 

21 

3.9E-11 

5.7E-07 

2.2E-f01 


17 

4.7E-02 

3.7E-02 

22 

3.5E-11 

4.8E-06 

1.2E-f01 

U.Ui 


33 

1.4E-02 

2.5E-02 

22 

5.8E-11 

3.3E-06 

6.1E-f00 


65 

2.8E-03 

l.lE-02 

23 

3.6E-11 

1.6E-06 

3.1E-f00 


129 

1.2E-03 

3.2E-03 

32 

6.5E-11 

3.7E-09 

1.6E-f00 


5 

1.6E-01 

4.0E-04 

21 

4.7E-11 

9.6E-17 

4.0E-f03 


9 

9.1E-02 

4.0E-04 

21 

3.1E-11 

5.6E-08 

2.2E-f03 

10-* 

17 

4.8E-02 

4.0E-04 

21 

6.5E-11 

l.OE-07 

1.2E-f03 

33 

2.5E-02 

4.0E-04 

22 

4.6E-11 

5.7E-07 

6.1E-f02 


65 

1.3E-02 

4.0E-04 

23 

4.2E-11 

1.3E-06 

3.1E-f02 


129 

6.3E-03 

4.0E-04 

30 

9.0E-11 

1.7E-06 

1.6E-f02 


5 

1.6E-01 

4.0E-06 

21 

4.7E-11 

7.2E-17 

4.0E-f05 


9 

9.1E-02 

4.0E-06 

21 

3.1E-11 

2.9E-12 

2.2E-f05 

10-6 

17 

4.8E-02 

4.0E-06 

21 

6.5E-11 

3.5E-12 

1.2E-f05 

33 

2.5E-02 

4.0E-06 

22 

4.6E-11 

1.4E-10 

6.1E-f04 


65 

1.3E-02 

4.0E-06 

23 

4.3E-11 

5.0E-10 

3.1E-f04 


129 

6.3E-03 

4.0E-06 

30 

9.0E-11 

3.3E-09 

1.6E-f04 


5 

1.6E-01 

4.0E-08 

21 

4.7E-11 

9.6E-17 

4.0E-f07 


9 

9.1E-02 

4.0E-08 

21 

3.1E-11 

l.lE-15 

2.2E-f07 

10-6 

17 

4.8E-02 

4.0E-08 

21 

6.5E-11 

4.3E-13 

1.2E-f07 

33 

2.5E-02 

4.0E-08 

22 

4.6E-11 

9.5E-12 

6.1E-f06 


65 

1.3E-02 

4.0E-08 

23 

4.3E-11 

1.3E-10 

3.1E-f06 


129 

6.3E-03 4.0E-08 30 9.0E-11 5.8E-12 

36 of 52 

American Institute of Aeronautics and Astronautics 

1.6E-f06 



Table 4. Constant Wave Speed - Shock - Even Number of Cells 


100 . 


1.00 


0.01 


10 -* 


10-6 


10 "® 


cells 

\u- Ue\ 

I'^Roe '^e\ 

iter. 

L 2 

^2, Roe 

<1 

6 

6.9E-08 

6.5E-07 

4 

2.9E-14 

4.2E-13 

3.3E-03 

12 

1.9E-08 

1.4E-07 

4 

6.3E-13 

1.6E-12 

1.7E-03 

24 

4.8E-09 

2.9E-08 

4 

4.1E-12 

6.6E-12 

8.3E-04 

48 

1.2E-09 

6.2E-09 

4 

4.4E-11 

5.2E-11 

4.2E-04 

96 

3.0E-10 

1.4E-09 

5 

5.0E-12 

2.6E-12 

2.1E-04 

192 

7.3E-11 

3.2E-10 

5 

3.8E-11 

5.5E-11 

l.OE-04 

6 

7.8E-04 

6.5E-03 

8 

6.1E-12 

6.4E-10 

3.3E-01 

12 

2.1E-04 

1.5E-03 

9 

3.1E-11 

l.lE-09 

1.7E-01 

24 

5.3E-05 

3.0E-04 

11 

7.5E-11 

3.1E-10 

8.3E-02 

48 

1.3E-05 

6.2E-05 

15 

5.4E-11 

8.7E-11 

4.2E-02 

96 

3.3E-06 

1.4E-05 

22 

3.9E-11 

4.5E-11 

2.1E-02 

192 

8.3E-07 

3.2E-06 

34 

5.1E-11 

5.4E-11 

l.OE-02 

6 

3.8E-08 

1.4E-02 

12 

4.6E-11 

1.8E-04 

3.3E+01 

12 

2.0E-08 

1.4E-02 

13 

2.6E-11 

1.3E-03 

1.7E+01 

24 

5.7E-05 

1.4E-02 

16 

2.5E-11 

2.5E-04 

8.3E+00 

48 

1.3E-03 

1.2E-02 

20 

5.0E-11 

8.5E-06 

4.2E+00 

96 

1.5E-03 

5.5E-03 

28 

2.7E-11 

5.0E-08 

2.1E+00 

192 

4.4E-04 

l.lE-03 

40 

8.9E-11 

1.7E-11 

l.OE+00 

6 

3.8E-08 

1.6E-04 

12 

4.5E-11 

1.2E-04 

3.3E+03 

12 

1.9E-08 

l.lE-04 

13 

2.1E-11 

2.7E-04 

1.7E+03 

24 

9.4E-09 

1.2E-04 

16 

2.9E-11 

2.4E-04 

8.3E+02 

48 

4.7E-09 

1.3E-04 

20 

5.5E-11 

3.0E-05 

4.2E+02 

96 

2.3E-09 

1.3E-04 

27 

3.0E-11 

2.1E-06 

2.1E+02 

192 

1.2E-09 

1.3E-04 

37 

6.4E-11 

7.6E-06 

l.OE+02 

6 

3.8E-08 

2.0E-06 

12 

4.5E-11 

3.5E-06 

3.3E+05 

12 

1.9E-08 

4.2E-07 

13 

2.1E-11 

8.9E-06 

1.7E+05 

24 

9.4E-09 

8.7E-07 

16 

2.9E-11 

9.0E-06 

8.3E+04 

48 

4.7E-09 

1.3E-06 

20 

5.5E-11 

5.8E-07 

4.2E+04 

96 

2.3E-09 

1.3E-06 

27 

3.0E-11 

8.5E-07 

2.1E+04 

192 

1.2E-09 

1.3E-06 

37 

6.4E-11 

1.5E-08 

l.OE+04 

6 

3.8E-08 

5.4E-08 

12 

4.5E-11 

4.4E-07 

3.3E+07 

12 

1.9E-08 

2.1E-07 

13 

2.1E-11 

2.2E-06 

1.7E+07 

24 

9.4E-09 

6.8E-09 

16 

2.9E-11 

3.9E-07 

8.3E+06 

48 

4.7E-09 

l.lE-08 

20 

5.5E-11 

l.OE-07 

4.2E+06 

96 

2.3E-09 

1.2E-08 

27 

3.0E-11 

7.7E-08 

2.1E+06 

192 

1.2E-09 

1.3E-08 

37 

6.4E-11 

9.9E-09 

l.OE+06 
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Table 5. 

V 

End Point Wave Speed - 

cells \u-Ue\ \uRoe-Ue 

Shock - 

iter. 

Even Number of Cells - (j> = 1 

L 2 L2,Roe Ax/v 


6 

6.9E-08 

6.5E-07 

4 

3.3E-14 

4.2E-13 

3.3E-03 


12 

1.9E-08 

1.4E-07 

4 

4.8E-13 

1.6E-12 

1.7E-03 

100. 

24 

4.8E-09 

2.9E-08 

4 

3.8E-12 

6.6E-12 

8.3E-04 

48 

1.2E-09 

6.2E-09 

4 

4.4E-11 

5.2E-11 

4.2E-04 


96 

3.0E-10 

1.4E-09 

5 

1.3E-12 

2.6E-12 

2.1E-04 


192 

7.5E-11 

3.2E-10 

5 

4.5E-11 

5.5E-11 

l.OE-04 


6 

8.0E-04 

6.5E-03 

8 

6.4E-12 

6.4E-10 

3.3E-01 


12 

2.1E-04 

1.5E-03 

9 

3.1E-11 

l.lE-09 

1.7E-01 


24 

5.3E-05 

3.0E-04 

11 

7.5E-11 

3.1E-10 

8.3E-02 

i.OO 


48 

1.3E-05 

6.2E-05 

15 

5.4E-11 

8.7E-11 

4.2E-02 


96 

3.3E-06 

1.4E-05 

22 

3.9E-11 

4.5E-11 

2.1E-02 


192 

8.3E-07 

3.2E-06 

34 

5.1E-11 

5.4E-11 

l.OE-02 


6 

l.OE-01 

1.4E-02 

11 

9.6E-11 

2.8E-04 

3.3E+01 


12 

5.2E-02 

1.4E-02 

12 

9.6E-11 

1.5E-03 

1.7E+01 


24 

2.6E-02 

1.4E-02 

16 

1.9E-11 

2.5E-04 

8.3E+00 

U.Ui 


48 

9.0E-03 

1.2E-02 

20 

7.8E-11 

8.5E-06 

4.2E+00 


96 

2.3E-03 

5.5E-03 

28 

3.1E-11 

5.0E-08 

2.1E+00 


192 

5.0E-04 

l.lE-03 

40 

9.0E-11 

1.7E-11 

l.OE+00 


6 

l.OE-01 

l.OE-04 

11 

9.6E-11 

1.6E-04 

3.3E+03 


12 

5.2E-02 

1.7E-04 

12 

9.2E-11 

3.5E-04 

1.7E+03 

10-* 

24 

2.6E-02 

1.2E-04 

16 

2.0E-11 

2.4E-04 

8.3E+02 

48 

1.3E-02 

1.3E-04 

20 

4.6E-11 

3.0E-05 

4.2E+02 


96 

6.4E-03 

1.3E-04 

27 

2.8E-11 

2.1E-06 

2.1E+02 


192 

3.2E-03 

1.3E-04 

37 

6.3E-11 

7.6E-06 

l.OE+02 


6 

l.OE-01 

4.4E-07 

11 

9.6E-11 

3.8E-06 

3.3E+05 


12 

5.2E-02 

2.5E-06 

12 

9.2E-11 

1.5E-05 

1.7E+05 

10-6 

24 

2.6E-02 

8.7E-07 

16 

2.0E-11 

9.0E-06 

8.3E+04 

48 

1.3E-02 

1.3E-06 

20 

4.6E-11 

5.8E-07 

4.2E+04 


96 

6.4E-03 

1.3E-06 

27 

2.8E-11 

8.5E-07 

2.1E+04 


192 

3.2E-03 

1.3E-06 

37 

6.3E-11 

1.5E-08 

l.OE+04 


6 

l.OE-01 

2.7E-07 

11 

9.6E-11 

1.3E-06 

3.3E+07 


12 

5.2E-02 

3.1E-07 

12 

9.2E-11 

8.2E-06 

1.7E+07 

10"® 

24 

2.6E-02 

6.8E-09 

16 

2.0E-11 

3.9E-07 

8.3E+06 

48 

1.3E-02 

l.lE-08 

20 

4.6E-11 

l.OE-07 

4.2E+06 


96 

6.4E-03 

1.2E-08 

27 

2.8E-11 

7.7E-08 

2.1E+06 


192 

3.2E-03 1.3E-08 37 6.3E-11 9.9E-09 

38 of 52 

American Institute of Aeronautics and Astronautics 

l.OE+06 



Table 6. 

V 

End Point Wave Speed - Shock - 

cells \u — Ue\ \uRoe — Ue\ Iter. 

Even Nnmber of Cells - <j> = 4’{4’) 
L 2 L2,Roe l^x/v 


6 

6.6E-08 

6.5E-07 

4 

4.4E-14 

4.2E-13 

3.3E-03 


12 

1.9E-08 

1.4E-07 

4 

4.0E-13 

1.6E-12 

1.7E-03 

100 

24 

4.8E-09 

2.9E-08 

4 

4.2E-12 

6.6E-12 

8.3E-04 

48 

1.2E-09 

6.2E-09 

4 

4.5E-11 

5.2E-11 

4.2E-04 


96 

3.0E-10 

1.4E-09 

5 

3.0E-12 

2.6E-12 

2.1E-04 


192 

7.6E-11 

3.2E-10 

5 

5.3E-11 

5.5E-11 

l.OE-04 


6 

7.6E-04 

6.5E-03 

8 

6.2E-12 

6.4E-10 

3.3E-01 


12 

2.1E-04 

1.5E-03 

9 

3.1E-11 

l.lE-09 

1.7E-01 


24 

5.3E-05 

3.0E-04 

11 

7.5E-11 

3.1E-10 

8.3E-02 

i.OO 


48 

1.3E-05 

6.2E-05 

15 

5.4E-11 

8.7E-11 

4.2E-02 


96 

3.3E-06 

1.4E-05 

22 

3.9E-11 

4.5E-11 

2.1E-02 


192 

8.3E-07 

3.2E-06 

34 

5.1E-11 

5.4E-11 

l.OE-02 


6 

9.3E-02 

1.4E-02 

11 

3.0E-11 

2.8E-04 

3.3E+01 


12 

4.7E-02 

1.4E-02 

12 

5.5E-11 

1.5E-03 

1.7E+01 


24 

2.3E-02 

1.4E-02 

16 

3.8E-11 

2.5E-04 

8.3E+00 

U.Ui 


48 

6.1E-03 

1.2E-02 

20 

6.1E-11 

8.5E-06 

4.2E+00 


96 

1.7E-03 

5.5E-03 

28 

2.8E-11 

5.0E-08 

2.1E+00 


192 

4.6E-04 

l.lE-03 

40 

9.0E-11 

1.7E-11 

l.OE+00 


6 

9.3E-02 

l.OE-04 

11 

3.0E-11 

1.6E-04 

3.3E+03 


12 

4.7E-02 

1.7E-04 

12 

4.0E-11 

3.5E-04 

1.7E+03 

10-* 

24 

2.3E-02 

1.2E-04 

16 

1.9E-11 

2.4E-04 

8.3E+02 

48 

1.2E-02 

1.3E-04 

20 

4.5E-11 

3.0E-05 

4.2E+02 


96 

5.8E-03 

1.3E-04 

27 

2.7E-11 

2.1E-06 

2.1E+02 


192 

2.9E-03 

1.3E-04 

37 

6.2E-11 

7.6E-06 

l.OE+02 


6 

9.3E-02 

4.4E-07 

11 

3.0E-11 

3.8E-06 

3.3E+05 


12 

4.7E-02 

2.5E-06 

12 

4.0E-11 

1.5E-05 

1.7E+05 

10-6 

24 

2.3E-02 

8.7E-07 

16 

1.9E-11 

9.0E-06 

8.3E+04 

48 

1.2E-02 

1.3E-06 

20 

4.5E-11 

5.8E-07 

4.2E+04 


96 

5.8E-03 

1.3E-06 

27 

2.7E-11 

8.5E-07 

2.1E+04 


192 

2.9E-03 

1.3E-06 

37 

6.2E-11 

1.5E-08 

l.OE+04 


6 

9.3E-02 

2.7E-07 

11 

3.0E-11 

1.3E-06 

3.3E+07 


12 

4.7E-02 

3.1E-07 

12 

4.0E-11 

8.2E-06 

1.7E+07 

10-6 

24 

2.3E-02 

6.8E-09 

16 

1.9E-11 

3.9E-07 

8.3E+06 

48 

1.2E-02 

l.lE-08 

20 

4.5E-11 

l.OE-07 

4.2E+06 


96 

5.8E-03 

1.2E-08 

27 

2.7E-11 

7.7E-08 

2.1E+06 


192 

2.9E-03 1.3E-08 37 6.2E-11 9.9E-09 
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Table 7. Constant Wave Speed - Boundary Layer - Odd Number of Cells 


V 


100 . 


1.00 


0.01 


10-4 


10 -® 


10 "® 


cells 

U - Me 

I'URoe 'Ue\ 

5 

l.OE-07 

2.3E-07 

9 

3.3E-08 

6.0E-08 

17 

9.5E-09 

1.9E-08 

33 

2.5E-09 

6.8E-09 

65 

6.6E-10 

2.1E-09 

129 

1.7E-10 

6.0E-10 

5 

8.6E-04 

1.6E-03 

9 

2.9E-04 

6.3E-04 

17 

8.5E-05 

2.1E-04 

33 

2.3E-05 

7.0E-05 

65 

5.9E-06 

2.1E-05 

129 

1.5E-06 

5.9E-06 

5 

2.2E-02 

8.7E-03 

9 

1.8E-02 

1.3E-02 

17 

1.2E-02 

9.6E-03 

33 

6.2E-03 

5.2E-03 

65 

2.5E-03 

3.4E-03 

129 

8.3E-04 

3.0E-03 

5 

l.OE-03 

l.OE-01 

9 

9.9E-04 

3.5E-04 

17 

9.1E-04 

4.8E-04 

33 

8.0E-04 

3.3E-04 

65 

6.7E-04 

2.8E-04 

129 

5.5E-04 

3.6E-04 

5 

2.7E-04 

2.8E-02 

9 

1.5E-04 

1.9E+10 

17 

8.0E-05 

1.2E+03 

33 

4.4E-05 

3.4E-06 

65 

2.7E-05 

3.4E-06 

129 

1.4E-05 

2.5E-06 

5 

2.7E-04 

2.7E-02 

9 

1.5E-04 

3.8E+17 

17 

7.9E-05 

1.4E+05 

33 

4.1E-05 

1.3E+20 

65 

2.1E-05 

3.6E+74 

129 

l.OE-05 

1.6+141 


iter. 

L 2 

^2, Roe 

3 

8.4E-11 

2.5E-10 

4 

2.0E-13 

l.OE-12 

4 

1.4E-12 

3.0E-12 

4 

1.2E-11 

1.6E-11 

5 

7.5E-13 

4.7E-13 

5 

l.OE-11 

9.4E-12 

7 

3.5E-11 

7.2E-10 

8 

3.5E-11 

1.8E-09 

10 

1.4E-11 

1.3E-10 

12 

6.8E-11 

1.5E-10 

17 

2.8E-11 

3.6E-11 

21 

4.6E-14 

2.2E-11 

26 

3.1E-11 

4.4E-08 

27 

5.3E-11 

2.0E-07 

28 

l.OE-11 

1.3E-06 

27 

1.5E-11 

8.0E-06 

26 

6.9E-11 

2.9E-05 

26 

3.2E-11 

2.4E-06 

31 

1.5E-11 

l.OE-02 

31 

3.1E-11 

3.9E-04 

31 

7.9E-11 

1.3E-03 

32 

9.2E-11 

1.4E-04 

33 

5.7E-11 

2.1E-03 

34 

2.6E-11 

3.2E-03 

47 

6.8E-11 

4.3E-04 

48 

7.3E-11 

4.8E+10 

49 

6.9E-11 

3.6E+03 

50 

6.4E-11 

1.3E-08 

51 

6.1E-11 

1.7E-06 

38 

2.5E-11 

l.OE-04 

47 

7.0E-11 

4.8E-04 

48 

7.4E-11 

9.9E+17 

49 

7.0E-11 

6.4E+05 

50 

6.5E-11 

7.2E+20 

51 

6.1E-11 

3.9E+75 

52 

5.8E-11 

5.2+142 


40 of 52 


^xjv 

4.0E-03 

2.2E-03 

1.2E-03 

6.1E-04 

3.1E-04 

1.6E-04 

4.0E-01 

2.2E-01 

1.2E-01 

6.1E-02 

3.1E-02 

1.6E-02 

4.0E+01 

2.2E+01 

1.2E+01 

6.1E+00 

3.1E+00 

1.6E+00 

4.0E+03 

2.2E+03 

1.2E+03 

6.1E+02 

3.1E+02 

1.6E+02 

4.0E+05 

2.2E+05 

1.2E+05 

6.1E+04 

3.1E+04 

1.6E+04 

4.0E+07 

2.2E+07 

1.2E+07 

6.1E+06 

3.1E+06 

1.6E+06 
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Table 8. 

V 

End Point Wave Speed - Boundary Layer - 

cells \u — Ue\ \uRoe — Ue\ Iter. L 2 

Odd Number of Cells 

L2,Roe ^x/v 


5 

l.OE-07 

2.3E-07 

3 

8.4E-11 

2.5E-10 

4.0E-03 


9 

3.3E-08 

6.0E-08 

4 

9.6E-14 

l.OE-12 

2.2E-03 

100. 

17 

9.5E-09 

1.9E-08 

4 

9.8E-13 

3.0E-12 

1.2E-03 

33 

2.5E-09 

6.8E-09 

4 

1.2E-11 

1.6E-11 

6.1E-04 


65 

6.6E-10 

2.1E-09 

5 

l.lE-12 

4.7E-13 

3.1E-04 


129 

1.7E-10 

6.0E-10 

5 

6.2E-12 

9.4E-12 

1.6E-04 


5 

8.2E-04 

1.6E-03 

7 

3.9E-11 

7.2E-10 

4.0E-01 


9 

2.9E-04 

6.3E-04 

8 

3.5E-11 

1.8E-09 

2.2E-01 


17 

8.5E-05 

2.1E-04 

10 

1.4E-11 

1.3E-10 

1.2E-01 

i.OO 


33 

2.3E-05 

7.0E-05 

12 

6.8E-11 

1.5E-10 

6.1E-02 


65 

5.9E-06 

2.1E-05 

17 

2.8E-11 

3.6E-11 

3.1E-02 


129 

1.5E-06 

5.9E-06 

21 

8.3E-14 

2.2E-11 

1.6E-02 


5 

8.9E-02 

8.7E-03 

39 

9.5E-11 

1.7E-09 

4.0E+01 


9 

4.0E-02 

1.3E-02 

40 

5.9E-11 

1.6E-10 

2.2E+01 


17 

2.4E-03 

9.6E-03 

43 

6.5E-11 

1.2E-12 

1.2E+01 

O.Oi 


33 

3.1E-03 

5.2E-03 

33 

2.3E-11 

2.0E-08 

6.1E+00 


65 

1.7E-03 

3.4E-03 

27 

1.8E-11 

4.0E-06 

3.1E+00 


129 

8.0E-04 

3.0E-03 

26 

6.2E-11 

2.4E-06 

1.6E+00 


5 

1.7E-01 

4.1E-01 

49 

7.3E-11 

3.8E-02 

4.0E+03 


9 

7.8E-02 

2.8E-04 

49 

8.7E-11 

2.5E-06 

2.2E+03 

10-4 

17 

4.1E-02 

3.2E-04 

48 

7.3E-11 

1.7E-05 

1.2E+03 

33 

2.1E-02 

3.3E-04 

49 

8.5E-11 

1.7E-06 

6.1E+02 


65 

l.OE-02 

3.2E-04 

49 

5.2E-11 

2.6E-05 

3.1E+02 


129 

5.1E-03 

3.1E-04 

49 

8.7E-11 

3.4E-05 

1.6E+02 


5 

1.7E-01 

2.9E-02 

49 

7.3E-11 

3.3E-04 

4.0E+05 


9 

7.9E-02 

5.2E+11 

51 

8.0E-11 

1.6E+12 

2.2E+05 

10-® 

17 

4.2E-02 

8.6E+03 

52 

5.6E-11 

2.8E+04 

1.2E+05 

33 

2.2E-02 

3.4E-06 

51 

9.3E-11 

9.9E-09 

6.1E+04 


65 

l.lE-02 

3.5E-06 

52 

9.7E-11 

1.3E-06 

3.1E+04 


129 

5.5E-03 

3.5E-06 

53 

9.0E-11 

2.6E-06 

1.6E+04 


5 

1.7E-01 

2.7E-02 

49 

7.3E-11 

3.7E-04 

4.0E+07 


9 

7.9E-02 

7.9E+22 

51 

8.0E-11 

1.6E+23 

2.2E+07 

10"® 

17 

4.2E-02 

2.1E+13 

63 

9.0E-11 

8.4E+13 

1.2E+07 

33 

2.2E-02 

1.6E+27 

53 

6.3E-11 

1.3E+28 

6.1E+06 


65 

l.lE-02 

l.lE+78 

53 

9.7E-11 

1.5E+79 

3.1E+06 


129 

5.5E-03 4.8+152 54 8.7E-11 INF 
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Table 9. End Point Wave Speed - Boundary Layer - Odd Number of Cells - <j> = 4’i4') 


V 


100 


1.00 


0.01 


10-4 


10 -® 


10 "® 


cells 

|U - Mel 

\^Roe ^e| 

iter. 

^2 

^2, Roe 

Ax/m 

5 

9.7E-08 

2.3E-07 

3 

8.5E-13 

4.4E-11 

4.0E-03 

9 

3.3E-08 

6.0E-08 

3 

3.1E-12 

3.2E-10 

2.2E-03 

17 

9.5E-09 

1.9E-08 

3 

1.4E-11 

3.1E-10 

1.2E-03 

33 

2.5E-09 

6.8E-09 

3 

7.4E-11 

3.0E-10 

6.1E-04 

65 

6.6E-10 

2.1E-09 

4 

6.4E-13 

2.6E-13 

3.1E-04 

129 

1.7E-10 

6.0E-10 

4 

6.1E-12 

1.7E-12 

1.6E-04 

5 

7.9E-04 

1.6E-03 

5 

1.6E-11 

1.5E-08 

4.0E-01 

9 

2.9E-04 

6.3E-04 

5 

7.5E-11 

3.9E-07 

2.2E-01 

17 

8.5E-05 

2.1E-04 

6 

7.9E-12 

l.lE-08 

1.2E-01 

33 

2.3E-05 

7.0E-05 

7 

3.9E-12 

2.0E-10 

6.1E-02 

65 

5.9E-06 

2.1E-05 

8 

8.2E-12 

3.3E-11 

3.1E-02 

129 

1.5E-06 

5.9E-06 

9 

6.2E-11 

9.4E-11 

1.6E-02 

5 

1.3E-02 

8.5E-03 

21 

4.7E-11 

3.5E-04 

4.0E-f01 

9 

1.5E-02 

1.3E-02 

23 

5.6E-11 

1.7E-06 

2.2E-f01 

17 

1.6E-02 

9.6E-03 

58 

7. IE- 11 

4.2E-16 

1.2E-f01 

33 

6.7E-04 

5.2E-03 

24 

1.6E-11 

2.5E-09 

6.1E-H00 

65 

1.8E-03 

3.4E-03 

24 

3.1E-11 

1.8E-07 

3.1E-H00 

129 

8.0E-04 

3.0E-03 

24 

6.1E-12 

4.6E-07 

1.6E-f00 

5 

4.9E-02 

5.7E-02 

26 

6.5E-11 

6.5E-03 

4.0E-H03 

9 

2.7E-02 

2.7E-04 

33 

8.3E-11 

1.3E-05 

2.2E-f03 

17 

1.4E-02 

3.1E-04 

36 

5.0E-11 

1.2E-04 

1.2E-f03 

33 

7.1E-03 

3.4E-04 

37 

5.4E-11 

1.4E-04 

6.1E-H02 

65 

3.2E-03 

3.2E-04 

38 

7.8E-11 

1.2E-04 

3.1E-H02 

129 

1.6E-03 

3.1E-04 

36 

6.0E-11 

5.0E-05 

1.6E-f02 

5 

4.9E-02 

3.4E-02 

26 

6.5E-11 

1.7E-03 

4.0E-H05 

9 

2.8E-02 

2.2E+00 

37 

3.0E-11 

5.8E-f00 

2.2E-f05 

17 

1.5E-02 

5.8E+06 

40 

5.3E-11 

1.7E-f07 

1.2E-f05 

33 

7.5E-03 

3.4E-06 

68 

7.9E-11 

2.8E-09 

6.1E-H04 

65 

3.8E-03 

3.4E-06 

40 

9.2E-11 

4.0E-06 

3.1E-H04 

129 

1.9E-03 

3.3E-06 

41 

8.5E-11 

1.2E-05 

1.6E-f04 

5 

4.9E-02 

3.5E-02 

26 

6.5E-11 

1.7E-03 

4.0E-H07 

9 

2.8E-02 

2.0E+01 

36 

7.0E-11 

8.0E-f01 

2.2E-f07 

17 

1.5E-02 

6.3E+07 

42 

6.7E-11 

2.5E-f08 

1.2E-f07 

33 

7.5E-03 

2.5E+32 

43 

8.3E-11 

1.8E-f33 

6.1E-H06 

65 

3.8E-03 

2.9E-08 

44 

8.7E-11 

2.6E-07 

3.1E-H06 

129 

1.9E-03 

1.2E+60 

45 

9.4E-11 

3.7E-f61 

1.6E-f06 
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Table 10. Convergence Times for Qusai, One-Dimensional Nozzle 


Case 

Cells 

Total Time 

Recon. Time 

Iter. 

Recon. Time 
Iter. 

Recon. Time 
(Iter.)(Cell) 


36 

1.67 

1.05 

273 

0.003846 

0.0001068 

EPWS 

70 

1.30 

0.93 

119 

0.0078151 

0.0001116 

Exit Area = 20 

138 

2.63 

1.73 

117 

0.014786 

0.0001071 

Pstag = 10® 

274 

5.40 

3.85 

129 

0.029844 

0.0001089 


546 

13.79 

9.92 

164 

0.060487 

0.0001108 


36 

1.21 

0.49 

406 

0.001206 

0.0000335 

Roe 

70 

1.23 

0.50 

247 

0.002024 

0.0000289 

Exit Area = 20 

138 

2.63 

1.01 

268 

0.003768 

0.0000273 

Pstag = 10® 

274 

7.65 

2.78 

401 

0.006932 

0.0000253 


546 

23.65 

8.82 

627 

0.014067 

0.0000258 


VII. Appendix A — One-Dimensional Reconstruction with Quadrature 

The methods in this section have not been fully explored in numerical tests but are a logical 
extension of the basic formulations in Sec. III. The methods require quadrature and tend not to 
be robust (or require careful application of limiters) in the limit of ^ 0. They are presented here 
to preserve ideas for future development of higher order methods. 


A. Multi-segment Wave Speed Approximation 

The multi-segment wave speed approximation breaks the interval into multiple segments. 

The wave speed across each segment is treated as a constant. This approach is essentially a solution 
by quadrature of Eq. 5. The linearized equation on segment N is 


dq d‘^q 

+ sn- 

dx dx^ 


(120) 


Following the derivation of Sec. A the transformed variable rj is defined. 


r]{x) = 


dr) 

dx 

d^7} 



( 121 ) 


The point x lies in the segment [xv-i-i, xv] (segment A -|- 1). Segments are numbered consecutively 
from the interface xq moving to the endpoints xl or xr. The analytic solution at each segment 
interface is Cl continuous. In the limit of Aq ^ 0 the mapping defined by Eq. 7 becomes 


T] 

dr) 

dx 

a^T) 

dx‘^ 


X - XV + 
1 
0 


(122) 
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where r]i\f represents the value of the transformed variable through the first N segments in Eq. 121. 


N 


m = Yl 


n=l 


^n-1 

An— 1 


exp 


An— l(3^n 3?o) 


^n—l 


exp 


An— l(2^n— 1 2^o) 


t'n-l 


(123) 


Substitution of Eq. 121 (or Eq. 122 if Ajv = 0) into Eq. 120 yields the following. 


(fq 

dr]'^ 


sn f dx 


VN \dr] 


(124) 


The general solution to Eq. 124 is 


= W+(*)^(1-TO) 

= w + (s)„ (!)„(’)-’») 


(125) 


w + (s)„S(“P 


Ajv(x-xq) 

UN 


exp 


\n(xn—^o) 

UN 


The particular solution to Eq. 124 is derived as follows. Eirst, integrate both sides of Eq. 124, 
returning to independent variable x on the left hand side to simplify the integration. 


( dq'' 


_1JV fV 1 

fdx\ 

\dv, 

' particular 

l^N -^VN 

[dvJ 


= _^N fX XO 

un Jxn-xo ^ 


2\n{x') 

Un 




(126) 


_ Sn 
Ajv 


( 

Ajv(x— xo) 


(^exp 

VN 

— exp 


Xn{xn—xq) 

Un 


. N 


Note that 

particular, N 


= 0; consequently, the gradient is fully defined by the general 

solution for this case where s, A , and v are treated as a constants over the interval. Integrate both 
sides of Eq. 126 in the same manner to complete the derivation of the particular solution. 


^particular 


fV ( 

JrjN \dri J 


particular 


dvj' 


_ ^ rx-XQ 

~ XnJxn-xo\ ^ 


AAf(x') 

UN 


^ fx xo [ 1 _ 

Ajv Jxn—xq ' P 


exp 


Ajv(a;jv— jeq) 
Un 


S_N_ 

Ajv 


(x -xn)-^ (exp 


^n(x'-{xn-xo)) 
Un 

Xn(x-xn) 


^dx' 


Un 


dx' 
- 1 


(127) 




-Sexp 


Ajv(xjv— xo) 

( 

Ajv(x-xq) 


Ajv(xjv-Xq) 


Un 

(^exp 

Un 

— exp 

Un 

)\ 


SN 

\n 


(x — xtv) — exp 


Xn(xn-xq) 

Un 


(r?(x) - r?Aj) 


In the limit as Aq ^ 0 the particular solution in Eq. 127 goes to the expected result 


^particular 


SN 
2 UN 


(x - xn)" 


(128) 


The solution to Eq. 120 in the across the + 1 interval may now be expressed. 

Q {x^ — Qgeneral T Qparticular 

= qn + 

I Sn 
“T \ 




{x — xn) — exp 

Xn{xn-xo) 

Un 

(r?(x) - 77 aj) 


(129) 
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The requirement for Cl continuity across each segment yields the following relations in which the 
initial value of q and its gradient on segment N is equal to the respective terminal values from the 
previous segment. 


dq 

dx 


(In = 


N 


-.N-l, 


q" (xn) 
dq^-^{xN) 
dx 


(130) 

(131) 


The values for Aat and vm are algebraic functions of q^ (and possibly )• The reconstruction 

algorithm may now be implemented as follows. First, integrate toward xl- Initial values for go and 
(^) specified and corresponding values for Aq, t'O; and sq are computed. Segment 0 from xq 

to Xl is defined such that i?e(xi — xq) is of order 1 or smaller. Values for qi and are computed 

from Eq. 129 for N = 0. Again, updated values for Ai, and si are computed. The integration 
continues in this manner segment by segment until xl is reached or until |r/jv| exceeds a specified 
ceiling. Record the value of at this terminal location. Repeat the integration using the same 
initial values to xr and record the value of g|j for the terminal location to the right. Perturbations 
to initial values are introduced in a shooting algorithm until (g^ — qi)'^ + (g|j — gj?)^ < e. The initial 

values for gg and that satisfy the inequality are used to construct the interface flux. This 

algorithm enables a more accurate resolution of a steep transition between adjacent mesh points 
as compared to the constant wave speed approximation of Sec. III-A or the endpoint wave speed 
approximation of Sec. III-B. Shooting algorithms generally carry significant overhead in terms of 
convergence time and algorithm robustness. 


B. Linear Wave Speed Approximation 


The linear wave speed approximation is an interface centered formulation as is the constant wave 
speed approximation of Sec. III-A but it includes the gradient of wave speed to better approximate 
high gradient transitions through Aq = 0. Model Equation 5 is approximated over the interval 
[xl,xr] by 


(Ao + AiVA' 


dq d^q 

^ ~ ^0 XT 

dx dx^ 


(132) 


where A,,,o = (^)q. The source term s and diffusivity u are treated as constants across the 
interval [xr , xr] . The wave speed A is expanded as a linear function in x to better characterize the 
reconstruction of diffusion terms in the limit of large cell Reynolds numbers, Rcx-xq = A(x — xq) jv 
as will be highlighted subsequently. Linear (or higher) expansions of the source term and diffusivity 
may introduce better reconstruction at the expense of a more complex relaxation algorithm. Such 
possibilities have not yet been explored. A semi-analytic solution to Eq. 132 (some numerical 
quadrature is required in the case of non-zero source terms) is derived by introducing the following 
map from x to r/. 


T] 

dr) 

dx 

d?r) 

d^ 


rx—xo 

Jo 


exp 


(Aq + A 


x,0 2 



exp 


(Ao + A,,o^)^ 

Ao+Aa;,o(ai— xq) dr} 

UQ dx 


(133) 


Substitution of Eq. 133 into Eq. 132 yields 


d^g 

dif' 


so 

exp 

I'O 



Aq + Aa;^o 


X — Xo 


X — Xq 
KQ 


( 134 ) 


45 of 52 


American Institute of Aeronautics and Astronautics 



The general solution to Eq. 134 is 


Qgeneral 


qo + 
qo + 
qo + 



(Ao + A 


X 

x,0 2 



( 135 ) 


This formulation of a semi-analytic reconstruction requires quadrature to evaluate the transformed 
variable r/. The formulation of the particular solution is not developed in detail here but one can 
follow the example of previous sections in order to see that two such quadratures would be required. 
The overhead here is quite similar to that encountered with the shooting algorithm of the previous 
section. For the case with sq = 0 the reconstruction is implemented as detailed in Eqs. 15 - 20 
using the present formulation of r/. However, the Newton iteration in Eq. 19 must be set up as a 
function of two variables because r] is a function of both Ao(?o) and Ax,o((^') )• 


VIII. Appendix B — Vectors and Matrices for Mnltidimensional, 

Navier-Stokes Eqnations 


A. Perfect Gas 

The vector of conserved variables q and primitive variables q are defined: 


q 


q 


p 

pU 

pV 

pW 

pE 

P 

U 

V 

w 

T 


Transformations from primitive to conservative variables 
(T~^dq = dq) are defined 


T 


T-i 


1 

U 

V 

W 

E 


u 

p 

V 

p 

w 

p 




0 0 
P 0 

0 p 
0 0 
pU pV 

0 

1 

p 

0 

0 

u 

pCv 


(136) 


(137) 


(Tdq = dq) and conservative to primitive 


0 0 

0 0 

0 0 

P 0 

pW pC^ _ 

0 0 0 ' 

0 0 0 

A 0 0 

0 A 0 

V w 1 

pCv pCv pCy 


(138) 
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where a = \{U‘^ + V‘^ + VF^), E = e + a, and de = C^dT. 

The linearization employs matrices A and B defined below. 


dfi ^ 


0 

1 

0 

0 

0 

Pa — 

-pU + 2U 

-pv 

-pw 

p 

-uv 

V 

u 

0 

0 

-uw 

w 

0 

u 

0 

paU - UH 

-pu^ + R 

-puv 

-puw 

pu + u 


0 0 0 0 0 

„ 0 Ip 0 0 0 

ahi ~ 

T7^ = ^Bn = i^ 0 0 p 0 0 

0 0 0 p 0 

_ 0 \pU pV pW ^ 

Diagonalization of Ai = RiAiR^^ is implemented with 


2 

0 

0 

1 

1 

2U 

0 

0 

U + c 

U-c 

2V 

2c^ 

0 

V 

V 

2W 

0 

2c^ 

w 

w 

2d 

2Fc2 

2Wc^ 

H + cU 

H-cU 


C2- 

- pd 

pu 

pv 

pw 

-p 

- 

-V 

0 

1 

0 

0 

- 

-w 

0 

0 

1 

0 

Pd 

-Uc 

C-PU 

-pv 

-pw 

p 

Pd 

+ Uc 

-c-pU 

-pv 

-pw 

p 


Rr' = 


U 0 0 0 

0 [/ 0 0 

Ai = 0 0 [/ 0 

0 0 0 [7 + c 

0 0 0 0 

The corresponding matrices in the Y direction are: 


= Ao = 


0 

0 

0 

0 

U-c 


0 

0 

1 

0 

0 


-vu 

V 

u 

0 

0 


Pd-v^ 

-pu 

-pv + 2V 

-pw 

p 

(143) 

-vw 

0 

w 

V 

0 


pdV - VH 

-pvu 

-pv^ + ff 

-pvw 

pv + v 



d{§) 


0 0 0 0 0 

0 p 0 0 0 

= i^B22 = 0 0 |p 0 0 

0 0 0 p 0 

_ 0 pu yv pw ^ 
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Diagonalization of A 2 = R 2 A 2 R 2 ^ is implemented with 


R 2 = 


1 

2 I? 


R2' = 


Ao = 


The corresponding matrices in the Z direction are: 

0 0 0 


2 

0 


0 

1 



1 

2U 

0 


2c^ 

u 



U 

2V 

0 


0 


c 

V 

— c 

2W 

2c 

2 

0 

w 



W 

2d 21 Tc2 

2Uc^ 

H + cV 

H 

-cV 

— Pa 

pu 

pv 

pw 

-P 

-W 


0 


0 

1 


0 

-u 


1 


0 

0 


0 

- Fc 

- 

-pu 

c 

-pv 

-pw 

p 

+ Vc 

- 

-pu 

—c 

-pv 

-pw 

p 

V 

0 

0 

0 

0 




0 

V 

0 

0 

0 




0 

0 

V 

0 

0 




0 

0 

0 

F + 

c 0 




0 

0 

0 

0 

V - 

c 




( 145 ) 


(146) 


^f3_A _ 


-wu 

W 

0 

u 

0 

-wv 

0 

w 

V 

0 

Pa-W^ 

-pu 

-pv 

-PW + 2W 

p 

PaW - WH 

-pwu 

-pwv 

-pw^ + H 

pw + w 


(147) 


dti3 

<*(§) 


Rs = 


1 




0 

0 

0 

0 

0 




0 

P 

0 

0 

0 


’B 33 = 

= 

0 

0 

P 

0 

0 




0 

0 

0 

Ip 

0 




0 

pU 

pV 

IpW 

pCp 

Pr _ 


'3 ^ is 

implemented with 



2 


0 

0 


1 

1 

2U 

2c2 

0 


u 

u 

2V 


0 

2c2 


V 

V 

2W 


0 

0 


W + c 

w 

— c 

2d 

2Wc^ 

2Uc 

2 ff + cW 

H - 

cW 


(148) 


( 149 ) 


R3-1 = 


(3U PV PW -P 
-U 1 0 0 0 

-y 0 1 00 

Pa-Wc -PU -PV C-PW P 
Pa + Wc -pU -pV -c-pW P 


48 of 52 


American Institute of Aeronautics and Astronautics 



As 


VF 0 0 0 0 

0 W 0 0 0 

0 0 W 0 0 

0 0 0 W + c 0 

0 0 0 0 W-c 


For convenience, the product R ^BT is defined below. 


R_igT-i 



2 0 0 -/3 -f3 

0 2Pr 0 0 0 

0 0 2Pr 0 0 

-2 0 0 /3+|Pr ,3-|Pr 

-2 0 0 /3-|Pr ,3+|Pr 


(150) 


(151) 


B. Mixture of Thermally Perfect Gases 

The vector of conserved variables q and primitive variables q are defined: 


q 


q 


Ps 

pU 

pV 

pW 

pE 

pGy 

Ps 

U 

V 

w 

T 

T, 


(152) 


(153) 


Transformations from primitive to conservative variables (Tdq = dq) and conservative to prim- 
itive (T“^dq = dq) are defined 


T 


T-i 


^s,r 

u 

V 

w 

E 


1 

pCv,tr 


^s,r 

_U 

P 

_V_ 

P 

_1P 

P 

{a - etr) 


pCv 


0 0 0 0 0 

p 0 0 0 0 

0 p 0 0 0 

0 0 p 0 0 


pU pV pW pCy^tr pCv 


0 0 

0 

1 

p 

0 

0 

u 

pCv,tr 

0 


0 

0 

0 

1 

p 

0 

V 

pCv,tr 

0 


0 


0 

0 

0 

1 

p 

w 

pCv,tr 

0 


pCy^l] J 

0 

0 

0 

0 


1 

pCv,tr 

0 


0 

0 

0 

0 


1 

pCv,tr 

1 

pCv,V 


( 154 ) 
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The linearization employs matrices A and B defined below. 



U{5g,r - Cg) 

Cg 

0 

0 

0 

0 


Ir-U^ 

-(3U + 2U 

-^v 

-f3W 

/3 


II 

II 

. 

I 

-uv 

V 

u 

0 

0 

0 

dq 

-uw 

w 

0 

u 

0 

0 


JrU - UH 

-^7/2 + ff 

-i3UV 

-i3UW 

pu + u 



-Uev 

ev 

0 

0 

0 

u 


( 155 ) 


dhi 




r 

CsCr) 

0 

0 

0 

0 

0 




0 


tp 

0 

0 

0 

0 




0 


0 

p 

0 

0 

0 


V 


0 


0 

0 

p 

0 

0 

(156) 


^ 1 htr,r 

- - 

(1^)) 

IpU 

pV 

pW 

pcp 

PUr 

pcp 

PVv 



1 

■ - 

(fe)) 

0 

0 

0 

0 

1 



where = M / Mg if diffusion of mass is driven by gradients of mole fraction or = 1 if diffusion 
of mass is driven by gradients of mass fraction. The Schmidt number for species s is defined 
Scs = fi/ pDg. The Prandtl number for conduction of translational-rotational energy is defined 
Prt^ = pCp/fjtr and Prandtl number for conduction of vibrational-electronic energy is defined 
Pr^ = ixCp/fjy. The evaluation of the Jacobian for diffusion of energy requires a summation 
over the product of mass diffusion flux times species enthalpies which is represented by (^) = 
CsCshtr.s /Scg and by Cg^ghy, s /Scg. 

Diagonalization of Ai = RiAiR^^ is implemented with 



26 g^r 

0 

0 

Cs 

Cs 

0 


27/ 

0 

0 

U + C 

u-c 

0 

1 

2V 

2 c 2 

0 

V 

V 

0 


21T 

0 

2 c 2 

w 

w 

0 


4q, _ 2711 

P 

2VP : 

2WP 

H + cU 

H-cU 

20 

0 


0 

0 

0 

6v 

€y 

2 


!,r Cgjr 

^Ucg 

PVcg 

PWcg 

-PCg 

-pCg 


-V 

0 

1 

0 

0 

0 


-w 

0 

0 

1 

0 

0 

Jr — Uc 

C-PU 

-pv 

-pw 

P 

p 

Jr + Uc 

-C-PU 

-pv 

-pw 

P 

p 



PUey 

pVey 

PWCy 

— PCy C^ 

— pc, 


(157) 


7 / 0 0 0 0 0 

0 7 / 0 0 0 0 

0 0 7/ 0 0 0 

000 U+c 0 0 

0 0 0 0 U-c 0 

0 0 0 0 0 7/ 


( 158 ) 
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The corresponding matrices in the V direction are: 



■ V(S,,r - Cs) 

0 

Cs 

0 

0 

0 


-vu 

V 

u 

0 

0 

0 

df2 . 

-1- = -^2 = 

Pr-V^ 

-f3U 

-pv + 2V 

-pw 

p 

p 

dq 

-vw 

0 

IT 

V 

0 

0 


JrV - VH 

-pvu 

-/3T2 + H 

-pvw 

pv + v 

pv 


-Vcy 

0 

Cy 

0 

0 

V 



I^B22 = 



0 0 

P 0 
0 

0 0 
pU \pV 

0 0 


Diagonalization of A 2 = R 2 A 2 R 2 ^ is implemented with 


25s^r 0 0 Cs 

2U 0 2c^ U 

2V 0 0 V + c 

2W 2c2 0 W 


4a 2Wc^ 2Uc^ H + cV 

0 0 0 Ct, 

C^Ss,r - Cs% PUCs PVCs PWCs 


-IT 

0 

0 

1 

-u 

1 

0 

0 

Ir - Vc 

-pu 

c-pV 

-pw 

7r + Vc 

-pu 

-C-PV 

-pw 


PUCy 

PVCy 

PWCy 



0 0 
0 0 
0 0 
p 0 

fe 

0 0 


0 

0 

0 

0 


pCp^ 

Ptv 

pcp 

Pr^; 


Cs 0 

U 0 

V -c 0 

W 0 

ff-cV 

By 2 

-(3cs -^Cs 

0 0 

0 0 

(5 ^ 

(3 ^ 

— (3ey CP — (pSy 


A 2 


y 0 0 0 0 0 

0 y 0 0 0 0 

0 0 1/0 00 

000 v+c 0 0 

0 0 0 0 V-c 0 

0 0 0 0 0 V 


The corresponding matrices in the Z direction are: 


W{5s,r - Cs) 

0 

0 

Cs 

0 

0 

-WU 

IT 

0 

U 

0 

0 

-ITT 

0 

IT 

V 

0 

0 

Ir-W^ 

-pu 

-pv 

-PW + 2W 

p 

p 

prW - WH 

-pwu 

-pwv 

-PW^ + H 

pw + w 

pw 

-WCy 

0 

0 

€y 

0 

IT 


( 159 ) 


(160) 


(161) 


(162) 


( 163 ) 
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dh3 

i(§) 


= = 


{ds,r 0 


P 

0 

0 


0 

0 

p 

0 


0 

0 

0 

^p 


0 

0 

0 

0 


Rs = 


1 

2c2 


C^ 5 o r — 


0 

C-s^r 


0 

/3[/c, 


0 

PVcg 


€-v 

fiWcg 


€y 

-PCs 




htr,r / htr \'\ 

Sc^ \Sc /) 

^v,r / /l-u \ \ 

Scr- \ Sc / y 

pU pV 1 
0 0 

pW ^ 
0 0 

pCp 

Pr„ 

pCp 

Pr„ 

R3A3R3 ^ 

is implemented with 



■ 


0 0 

Cs 

Cs 

0 


2 U 

2c^ 0 

U 

U 

0 


2R 

0 2c2 

V 

V 

0 


2PR 

0 0 

PR + c 

PR-c 

0 

4a- ^ 

/3 

2PRc2 217C 

H + cW 

H-cW 

20 

p 


2 


-u 

1 

0 

0 

0 

0 

-V 

0 

1 

0 

0 

0 

7 ,. — PRc 

-pu 

-pv 

c-pW 

p 

p 

7 ^ + PRc 

-pu 

-pv 

-C-PW 

p 

p 

Cypr 

pUe^ 

pVe, 

PWe, 

-j3e^ 

(P — pe^ 


A3 


R" 0 0 0 0 0 

0 PR 0 0 0 0 

0 0 PR 0 0 0 

0 0 0 PR + c 0 0 

0 0 0 0 PR-c 0 

0 0 0 0 0 PR 


(164) 


(165) 


(166) 
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